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Abstract 


Multi-dimensional Upwind Fluctuation Splitting Scheme with Mesh 
Adaption for Hypersonic Viscous Flow 

By 

William Alfred Wood III 

Committee Chairman: Bernard Grossman 
Aerospace Engineering 

(ABSTRACT) 

A multi-dimensional upwind fluctuation splitting scheme is developed and imple- 
mented for two-dimensional and axisymmetric formulations of the Navier-Stokes equa- 
tions on unstructured meshes. Key features of the scheme are the compact stencil, 
full upwinding, and non-linear discretization which allow for second-order accuracy 
with enforced positivity. Throughout, the fluctuation splitting scheme is compared 
to a current state-of-the-art finite volume approach, a second-order, dual mesh up- 
wind flux difference splitting scheme (DMFDSFV), and is shown to produce more 
accurate results using fewer computer resources for a wide range of test cases. The 
scalar test cases include advected shear, circular advection, non-linear advection with 
coalescing shock and expansion fans, and advection-diffusion. For all scalar cases the 
fluctuation splitting scheme is more accurate, and the primary mechanism for the 
improved fluctuation splitting performance is shown to be the reduced production of 
artificial dissipation relative to DMFDSFV. The most significant scalar result is for 



combined advection-diffusion, where the present fluctuation splitting scheme is able to 
resolve the physical dissipation from the artificial dissipation on a much coarser mesh 
than DMFDSFV is able to, allowing order-of-magnitude reductions in solution time. 
Among the inviscid test cases the converging supersonic streams problem is notable in 
that the fluctuation splitting scheme exhibits superconvergent third-order spatial ac- 
curacy. For the inviscid cases of a supersonic diamond airfoil, supersonic slender cone, 
and incompressible circular bump the fluctuation splitting drag coefficient errors are 
typically half the DMFDSFV drag errors. However, for the incompressible inviscid 
sphere the fluctuation splitting drag error is larger than for DMFDSFV. A Blasius 
flat plate viscous validation case reveals a more accurate ^-velocity profile for fluctu- 
ation splitting, and the reduced artificial dissipation production is shown relative to 
DMFDSFV. Remarkably the fluctuation splitting scheme shows grid converged skin 
friction coefficients with only five points in the boundary layer for this case. A viscous 
Mach 17.6 (perfect gas) cylinder case demonstrates solution monotonicity and heat 
transfer capability with the fluctuation splitting scheme. While fluctuation splitting 
is recommended over DMFDSFV, the difference in performance between the schemes 
is not so great as to obsolete DMFDSFV. The second half of the dissertation develops 
a local, compact, anisotropic unstructured mesh adaption scheme in conjunction with 
the multi-dimensional upwind solver, exhibiting a characteristic alignment behavior 
for scalar problems. This alignment behavior stands in contrast to the curvature 
clustering nature of the local, anisotropic unstructured adaption strategy based upon 
a posteriori error estimation that is used for comparison. The characteristic align- 
ment is most pronounced for linear advection, with reduced improvement seen for 
the more complex non-linear advection and advection-diffusion cases. The adaption 
strategy is extended to the two-dimensional and axisymmetric Navier-Stokes equa- 
tions of motion through the concept of fluctuation minimization. The system test 
case for the adaption strategy is a sting mounted capsule at Mach-10 wind tunnel 
conditions, considered in both two-dimensional and axisymmetric configurations. For 
this complex flowfleld the adaption results are disappointing since feature alignment 
does not emerge from the local operations. Aggressive adaption is shown to result 
in a loss of robustness for the solver, particularly in the bow shock/stagnation point 
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interaction region. Reducing the adaption strength maintains solution robustness but 
fails to produce significant improvement in the surface heat transfer predictions. 
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Chapter 1 
Introduction 


1.1 Objective 

The present dissertation seeks to develop a next-generation numerical aerothermody- 
namic predictive capability. High-fidelity fiowfield solutions are sought with increased 
accuracy and reduced analysis time, especially with respect to configuration changes. 
This is an effort to push the state of the art for aerothermodynamic analysis capability 
closer to being a usable tool for vehicle designers. 

Specifically, the multi-dimensional upwinding concepts of Sidilkoverfl, 2, 3] as 
applied to the Euler equations of inviscid, perfect gas flows are incorporated into a 
consistent treatment of viscous and heating terms to solve the Navier-Stokes equa- 
tions of gas dynamics. A detailed analysis of the base algorithm is performed and 
applications to several validation cases are conducted. Also, aggressive unstructured 
mesh-adaption strategies are investigated in conjunction with the fluctuation splitting 
scheme. 


1.2 Motivation 

Current proposals for future United States access-to-space systems incorporate reus- 
able trans-atmospheric flight vehicles, performing a role similar to the space shuttle 
orbiter, seeking to minimize recurring costs. For efficient design and operation of these 
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vehicles aerothermodynamic performance predictions are required to high accuracy 
and in a timely fashion. 

The current computational capability to provide these aerothermodynamic per- 
formance predictions for complex vehicles is severely limited by total solution times, 
which include both domain discretization and fiowfield evolution, measured in months. 
In order for high-fidelity aerothermodynamic predictive tools to play an active role 
in the design phase a leap in responsiveness must be made over the current state of 
the art methods, represented by dimensionally-split approximate Riemann solvers. 

An effective numerical aerothermodynamic predictive tool must be able to provide 
rapid analysis of vehicle aerodynamics and control surface effectiveness. To do so 
fiowfield features, such as embedded shocks, shear layers, and boundary layers, must 
be accurately modeled and solution domains must be easily generated for complex 
shapes, allowing rapid re-calculation with respect to geometry changes. 

The emphasis in aerothermodynamics on heat transfer and high-speed flows pro- 
vides additional challenges with regard to the accurate resolution of boundary layers, 
loss of accuracy with highly-stretched meshes, and convergence slowdowns associated 
with a disparity in information speeds between nearly stagnated boundary layer flow 
and hypersonic shock layers. 

To make the next generation of computational aerothermodynamic predictive tools 
responsive to the design phase a truly multi-dimensional, robust Navier-Stokes solver 
based on general unstructured domains is required. Advanced, aggressive conver- 
gence acceleration methods and the exploitation of distributed or massively parallel 
architectures are mandatory. 


1.3 Background 

1.3.1 Fluid Dynamics Algorithms 


The twin goals of improving both accuracy and efficiency have long been objectives of 
the computational fluid dynamics developers. Accuracy is typically addressed through 
the spatial discretization, with the Van Leer ‘Ultimate’ series of papers[4, 5] being a 
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driving force throughout the 1970’s. In the domain of hypersonic applications, effi- 
ciency has been improved through the temporal discretization, with MacCormack’s 
work standing out on explicit predictor/corrector [6], implicit line Gauss-Seidel[7, 8], 
and approximate factorization schemes[9], along with a notable extension by Can- 
dler to the LU-SGS scheme[10]. Several industrial solvers, such as GASPfll], have 
emerged from this finite volume legacy, such as the successful CFL3D code[12]. 

In his landmark paper, Roe[13], building on the work of Godunov[14], introduced 
an upwind, approximate- Riemann-problem solution technique for the one-dimensional 
Euler equations. The Roe scheme has been extended on structured meshes to multiple 
dimensions by decomposing the domain into locally one-dimensional problems aligned 
with the computational coordinates. Viscous terms have been incorporated and in a 
decade’s time the structured Roe scheme 1 for Navier-Stokes equations set the standard 
for continuum aerothermodynamics[15, 16, 17, 18]. 

Unstructured schemes have since been developed, seeking both much greater flex- 
ibility for complex geometries and a reduction in preferential solution directions 
aligned with computational coordinates. Notable contributions for the unstructured 
Roe schemes have been made by Barth[19, 20. 21, 22] and Whitaker[23, 24], These 
approaches, however, still rely upon a locally one-dimensional, dimensionally-split, 
algorithm at cell faces. 

A different tack has been taken by Fey[25, 26], who has advocated the method of 
transport for a true multi-dimensional treatment of the Euler equations. This method 
is more an extension of flux vector splitting concepts, an approach known to be more 
dissipative than flux difference split,ting[17]. 

The advent of multi-dimensional linear advection schemes, termed fluctuation 
splitt,ing[27, 28, 29, 30, 31, 32], induced Roe to revisit his one-dimensional scheme. 
The fluctuation splitting formulation is based on an unstructured triangulation and 
is local to each cell, suitable for application on massively parallel computers. A 
two-dimensional analog to the Roe scheme for the Euler equations has been devel- 
oped based on characteristic wave decompositions[33, 34, 35, 36, 37, 38, 39, 40, 41]. 
Unfortunately, these wave decomposition schemes have been of limited robustness 

1 Includes extensions of the base Roe scheme to higher-order accurate discretizations. 
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and have not reached their full potential. Applications have also been made to the 
shallow- water equations [42], 

Taking the wave decomposition idea to the extreme, some authors have applied 
fluctuation splitting to Boltzmann schemes for the Euler equations[43, 44, 45, 46]. 
These gas kinetic schemes are very expensive to run and show very modest improve- 
ments over dimensionally-split finite volume[47]. 

Diffusion terms have been incorporated into scalar fluctuation splitting schemes[21, 
48] and viscous terms have been added to the wave-decomposition models to solve 
Navier-Stokes problems[49, 50, 51, 52]. While the Galerkin approach to viscous terms 
appears compatible with fluctuation splitting for convective terms, these schemes suf- 
fer from a lack of robustness. Caret,te[51] states, “However, this scheme is not sat- 
isfactory at present in terms of robustness and convergence, and improvements in 
this respect is still subject [sic] of current research.” Even for the classical flat plate 
boundary layer problem Tomaich[50] reports, “...the agreement with the Blasius 
solution is rather poor.” 

Sidilkover, who along with Roe established a strong link between fluctuation split- 
ting and upwind flux difference splitting finite volume[3], has proposed an alternative 
multi-dimensional treatment of the Euler system [1, 2], Rather than performing a 
wave decomposition to decouple the Euler equations, Sidilkover employs fluctuation 
splitting to treat the system of equations as a whole unit. In doing so he has been 
able to apply efficient multigrid [53, 54, 55] solution strategies to shock reflection and 
channel flow problems. Advantages of Sidilkover’s method for the Euler equations in- 
clude: arbitrary triangulations, stability of Gauss-Seidel relaxation on high resolution 
discretization, compact stencil, second-order accuracy, and rotationally invariant ar- 
tificial dissipation. Among these, Sidilkover[55] claims, “The fundamental advantage 
of this approach is that it leads to a scheme that combines high-resolution and good 
stability properties.” 

Sidilkover’s fluctuation splitting scheme for the Euler equations has the promise 
to satisfy many of the cutting-edge aerodynamic requirements. R provides a truly 
multi-dimensional treatment of the governing system of equations in an unstructured, 
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compact, high-order algorithm suitable to distributed and massively parallel compu- 
tation. Being stable with respect to Gauss-Seidel relaxation, as opposed to Runge- 
Kutta[56] marching, opens the possibility of greatly improved multigrid convergence 
rates. Robustness appears to be improved relative to the wave decomposition models. 

What is lacking from the scheme to be an aerothermodynamic tool are the exten- 
sion to the Navier-Stokes equations, rigorous evaluation of robustness and accuracy 
relative to finite volume schemes on complex geometries, analysis of the effect of grid 
stretching, and critical determination of the accuracy and efficiency for heat transfer 
and skin friction calculations. 


1.3.2 Mesh Adaption Strategies 

Solution adaptive remeshing techniques have been utilized with some success for hy- 
personic flows on structured domains. A simple, effective approach developed by 
Gnoffo[57] utilizes a spring analogy energy minimization to align the bow shock and 
cluster to the boundary layer. This approach works very well for entry forebodies, for 
which it was developed, but is more difficult to apply to complex vehicle shapes. The 
method also is unresponsive to embedded shocks or other shock-layer flow features. 

Harvey[58, 59, 60] has developed a mesh adaption technique that is sensitive to 
shock-layer features to obtain parabolized Navier-Stokes solutions over simple config- 
urations, e.g. cones, using a spring analogy based on Gnoffo’s work. Unfortunately, 
defining relative clustering strengths for various flow features proved difficult, and a 
damaging lack of robustness is shown for three-dimensional leeside fiowfields[59]. 

Mesh adaption on unstructured domains offers a significant benefit over structured 
mesh adaption — the ability to insert and delete nodes. Much of the research in this 
area has gone into global remeshing using isotropic cells[61, 62, 63, 64, 65, 66, 67]. 
Grids composed of (nearly) isotropic cells quickly become prohibitively large for hy- 
personic applications, where the capture of essentially one-dimensional flow features, 
such as shocks, results in refinement in all three dimensions. These methods typically 
employ gradient clustering, using a second derivative check on some or all of the de- 
pendent variables to define clustering strengths. This sort of clustering is intended 
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Figure 1.1: Pictorial of “ideal” mesh for shock discontinuity (left) and mesh resulting 
from gradient based clustering (right). 


to reduce the interpolation error in a piecewise-linear data representation[68, 69], but 
is not necessarily driven by the flow physics, and can lead to excessive clustering or 
conflicting requirements in certain regions, such as a bow shock or stagnation point. 
Figure 1.1 presents an illustrative pictorial based on the results of Ait-Ali-Yahia[70] 
et al, where 18 cells were driven into the bow shock by gradient based clustering. A 
shock is pictured on the left side of the figure with an ideal mesh for a three-point 
stencil. On the right side is a mesh typically produced by gradient-based cluster- 
ing, which gathers more points in the vicinity of the discontinuity without providing 
sharper resolution. 

More recently, impressive results using anisotropic elements have been reported 
by Habashi[71, 72, 73, 74, 75, 76] et al. With that approach, all grid adaptions 
are local operations, as opposed to global remeshings. Highly-stretched elements are 
obtained, achieved by equating the interpolation error along each edge, again using 
a spring analogy minimization. The clustering is still driven by second derivatives in 
the solution. 

Roe[32, 77] has applied the concept of performing global mesh adaption via lo- 
cal node movements to a scalar advection problem using fluctuation splitting. His 
analysis reveals that a characteristic mesh results, with far fewer points required than 
gradient based clustering would imply. His method is based on the minimization of an 
objective function formed by taking derivatives of the fluctuation splitting scheme. A 
differentiable high-resolution linear scheme, which cannot be monotonic[56], is chosen. 
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For a complex fiowfield or for systems of equations it may not be possible to achieve 
a perfectly-aligned characteristic mesh, in which case the non-monotonic property 
would most likely be detrimental or even fatal to the solution. It is not clear how 
to extend the differentiability requirement to a high-resolution, non-linear monotonic 
scheme. 

Another attempt at mesh alignment has been performed by Trepanier[78] et al. 
While their scheme is unresponsive to characteristic lines, a wave decomposition model 
is used with an inviscid cell-centered finite volume solver to produce shock alignment 
of unstructured isotropic cells. Less success was demonstrated for shear alignment. 
One drawback of the method is the need to explicitly detect and classify flow features, 
which could hinder the extension to three dimensions. Similarly, Parikh[79] et al. use 
a wave decomposition in conjunction with a cell-centered finite volume solver to drive 
edge alignment, but with automated general feature detection. Unfortunately, their 
approach was able to provide only extremely modest benefits, in part because there 
was no mechanism for providing cell stretching. 


1.4 Course of Action 

A systematic approach to developing new methods is pursued. Current leading-edge 
technology is embraced and developed into a complete gas dynamics solver, applicable 
to reentry vehicles across their speed range. Throughout, the scientific method is 
followed, where new schemes are critically evaluated against the current state of the 
art. 

One-dimensional scalar advection is considered first, where the common ancestry 
of the upwind finite volume and fluctuation splitting paradigms is presented. Exten- 
sion to the incorporation of diffusive terms is presented, followed by the treatment 
of systems, for both the Euler and Navier-Stokes equations. The chapter culminates 
with the current state of the art in dimensionally-split schemes. Fluctuation splitting 
provides an identical discretization in one dimension. 

Two-dimensional scalar advection on unstructured domains is addressed next. The 
methodology for applying the dimensionally-split concepts in multiple dimensions is 
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shown. The different approach of fluctuation splitting is highlighted, providing a truly 
multi-dimensional treatment. A critical comparison between the upwind finite volume 
and fluctuation splitting schemes for scalar equations is performed. The extension 
to diffusive problems is developed, with emphasis placed on consistent incorporation 
into the fluctuation splitting framework while still maintaining formal second-order 
accuracy. 

Novel concepts for mesh adaption are developed analytically for two-dimensional 
scalar problems, based on the physics of the solution rather than gradient-based 
clustering. Emphasis is placed on concepts that can be extended to systems. Demon- 
strations for simple basic problems and more challenging advection-diffusion problems 
are performed. 

Systems in two dimensions are formulated for the Navier-Stokes equations, rep- 
resenting the crux of the development of the fluctuation splitting concepts. The 
treatment follows the lead of Sidilkover in treating the system as a cohesive unit, 
in contrast to the wave-model simplifications of Roe and Deconinck. Specialization 
to axisymmetric equations of fluid motion provide useful applications and consistent 
treatment of source terms. 

The mesh adaption strategies developed for scalar equations are reformulated for 
the fluid dynamics systems. The resolution of flowfield features is addressed, along 
with robustness and convergence of the adaption process. An automatic refinement 
procedure is sought whereby the solution is converged uniformly with mesh density. 

A critical aerothermodynamic evaluation of the fluctuation splitting algorithm and 
mesh adaption strategies are performed using the primary case of an entry capsule 
at Mach-10 wind tunnel conditions, for which prior experimental and computational 
data exists. 



Chapter 2 

One-Dimensional Analysis 


The current state of the art for solving the compressible Navier-Stokes equations, 
namely upwind fiux-difference-split finite volume schemes, is developed on non-uni- 
form meshes in one spatial dimension. Upwind flux difference splitting, in particular 
the Roe scheme[13, 80], is considered the most accurate scheme for compressible 
Navier-Stokes applications[17, 18], primarily because of the low levels of artificial 
dissipation introduced through the matrix dissipation model. The particular finite 
volume scheme considered in this dissertation can be described as node-based, median- 
dual mesh upwind Roe flux difference splitting finite volume with limited multi-di- 
mensional reconstruction, and will be abbreviated as DMFDSFV, for dual mesh flux 
difference splitting finite volume, throughout. 

The sections in this chapter progress from scalar advection to diffusion and then 
to combined advection-diffusion. Treatment of scalar equations is followed by the 
Euler and Navier-Stokes systems of equations. The fluctuation splitting approach is 
developed in parallel with the DMFDSFV analysis, and is shown to result in identical 
discretizations in one dimension. 

This chapter has a threefold purpose: to introduce DMFDSFV and fluctuation 
splitting basics, to develop the state of the art for the locally one-dimensional approx- 
imate Riemann solver used in the finite volume algorithms, and to serve as a prelude 
for the multi-dimensional analyses, where fluctuation splitting offers new capability 
over the DMFDSFV extensions. 
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Figure 2.1: One-dimensional finite volume domain. 


2.1 Domain 

In one dimension the domain considered is a discretization of the :r-axis, either with 
uniform or non-uniform spacing between grid nodes, which are indexed by i. Depen- 
dent variables are stored at the nodes. 

In the node-based finite volume context a median-dual control volume is con- 
structed about each node by defining a cell face halfway between adjacent nodes. 
This convention is depicted in Figure 2.1, with the cell faces referred to as the 
points. In the illustrative case of Figure 2.1, the i + \ node is a boundary point, and 
the corresponding cell extends only from i + \ to % + 1. The generalized volume of 
the median-dual cell is formed as x i+ i — x t _i_ for interior nodes or Xi+\ — x i+ i for 
the boundary point shown in Figure 2.1. If the solution at the boundary nodes is 
specified, then there are two fewer interior finite volumes to solve than the number 
of nodes. If the boundary nodes are updated by the interior scheme, then there is 
a one-to-one correspondence between the nodes and control volumes, plus two addi- 
tional numerical boundary conditions, which may be implemented with ghost nodes 
(z + 2 in Figure 2.1) or specified boundary fluxes. If the ghost node is used, it can be 
located co-incident with the physical boundary node, and does not need to have an 
associated control volume. A ghost node co-located with the physical boundary node, 
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Figure 2.2: One-dimensional fluctuation splitting domain. 


without an associated control volume extends more easily to multiple dimensions with 
unstructured grids than would a formulation based upon physically locating the ghost 
node. 

Notice that for the nodal distribution depicted, which has non-uniform spac- 
ing, the cell centroids, denoted by x in Figure 2.1, do not coincide with the nodes. 
Barth [22] states that nodal storage in this case, referred to as mass lumping in a 
finite element context, only alters the time accuracy of finite volume schemes, and 
not the steady state solutions. In the present dissertation unsteady problems will 
employ uniform grids, though stationary problems are free to use non-uniform node 
distributions. 1 

Notice also the piecewise-linear representation of data in the finite volume context. 
Discontinuous jumps in the dependent data are allowed at cell faces. 

Figure 2.2 depicts the discretization of the domain for the fluctuation splitting 
approach. The data is now continuous and piecewise linear over elements defined 
by the nodes. The centroids of the fluctuation splitting elements are at the same 
locations as the faces of the finite volume cells. No special definition of a boundary 

lr The use of uniform meshes to avoid degradation of time accuracy due to mass lumping is done for 
convenience. It is possible to compute cell centroids and use that location in calculations involving 
cell distances. Conversely, a cell-centered, rather than node based, data storage structure could be 
adopted for the finite volume discretization. 
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cell is required. The length of an element is, 

I'i.i — 1 %i -\- 1 (2.1) 

There are one fewer elements than nodes, and each element is associated with two 
nodes. Correspondingly, each interior node is associated with two elements, while 
each boundary node is associated with only one element. 

The data structure for fluctuation splitting shown in Figure 2.2 resembles a cell- 
centered finite volume layout, with a key difference that cell-centered finite volume 
stores the solution at cell centers, as opposed to at the nodes for fluctuation splitting. 
Additionally, it is emphasized that the fluctuation splitting data is C-0 continuous 
whereas cell-centered finite volume would have discontinuous jumps in the data at 
the cell interfaces. 


2.2 Scalar Advection 


A hyperbolic conservation law takes the form, 

U, + V-F =0 (2.2) 

where U is the vector of conserved variables and F is the flux 2 of these variables. 


Following Godunov[14], Eqn. 2.2 can be evaluated in an integral sense, 


/ u, in = - I 

Jn Jn 


V • F dQ 


Jn Jn 

If the control volume, £2, is fixed in time, then, 


Jn 


U t dn =S n U f 


(2.3) 


(2.4) 


where the over-bar indicates a cell-average value. Using the divergence theorem the 
flux term can be evaluated as, 


/ V -FdQ = J)F-ndT 
Jn Jr 


(2.5) 


2 The flux function F is a function of the dependent variables U, F = F(U). For some of the 
scalar cases F will also be a function of the independent variables, F = F(u, x, y). In such cases the 
flux will be in a variables-separable form, F = \(x,y)u. 
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where n is the outward unit normal to the control volume boundary, V. 

The one-dimensional scalar advection problem is obtained from Eqn. 2.2 as, 

u,t + F x =0 (2-6) 

which can be written for a control cell as, 

Snu t = - / F x dQ F left face bright face (2.7) 

Jn 

2.2.1 Linear Advection 

Linear advection is obtained from Eqn. 2.6 by choosing F = Xu. The advection speed, 
A, is taken to be constant. 
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The first-order upwind is then seen to be the same as a central distribution plus an 
artificial dissipation term, 

$u = y 5 > v ( 2 - 14 ) 

Second-order upwind is constructed following the MUSCL concept of Van Leer[5], 
where a linear reconstruction is performed on each finite volume. The numerical flux 
of Eqn. 2.11 is modified to be, 



A-l- 1 A | A — | A | 


(2.15) 


where U|_ is the reconstructed conserved variable on the left side of the cell face and 
ur is the reconstructed variable on the right side of the cell face. Following Barth[22], 
a limited reconstruction is performed on each cell as, 


u fo 


Ui + i'i(Vu)i ■ f fa 


(2.16) 


The gradient is evaluated as a central difference, 

(*«)< = (f (2.17) 

The limiter function, ip, is employed to provide monotonicity of the solution, based 
upon positivity arguments. The limiter takes the form, 


where, 


, , P 

ip - 

\Q 


U{±_ 1 U,; . 

P = 5 • <1 = ( Vm ) 


r i±k 


The more restrictive of the ± choices is used for ip. Some popular limiters are pre- 
sented in appendix A. 

The discrete numerical flux (Eqn. 2.15) expands to, 


fi 


i+i 


A+|A| ( . £ij + 1 . ^ A— |A| f 

— l * + ^^2S~ diV ) + ~^2~ ( i+i ~ ^ 


- ) (2.18) 


2 S, 


i+1 


T Fi+i | A | ^ ^i,i+l ( Pi X ! T? I I \ I 1 1pi + 1 x IT? Ill \ 

- —A iu + — — I - - — 5; + i(F— |A|u) 


Si 


S, 


i+l 
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leading to, 


Ri — —6iF + d’u + i?2 


where the second-order correction is, 


Riu = ^(|^i-i(F+AW-|*(F-|AM 


( 2 . 20 ) 


The residual (Eqn. 2.19) can be rearranged as, 

Ri = SiF - l [/, 2 - 2 iPiF^ + - t^p.) Ft 

O L Oi - 1 \ Oj+1 Oi-l J 


+ 2 tpiFi+i — £ij+ Fi + 2 + d>2u 

<->2+1 

where the artificial dissipation is now, 


( 2 . 21 ) 


— $u + y _ — Ut- 2 + — (£j,j+i — ^i_i,t) Uj_i 

(o 1 , ,, ti + A i’i to o \ 

( ^2— 1,2 q T H.i+l „ I Ui \Z-i,i + 1 f-i— li) '^2+1 


, Vi+l 
''2,2+1 “2+2 
^2+1 


( 2 . 22 ) 


On a uniform grid and without limiting, the second-order residual (Eqn. 2.21) 
reduces to a low-truncation-error central difference minus fourth-order dissipation, 

^ = o (~T’j!-2 + 6E)_i — 6E) + i + Fj +2 ) + y (—Ui- 2 + — 6uj + 4iq + i — tq +2 ) 

o o 

(2.23) 

Fluctuation splitting 


In the fluctuation splitting framework Eqn. 2.7 is evaluated over each domain element, 
without recourse to the divergence theorem. The element fluctuation is defined as, 


Snu t = (Pe = ~ / F x (El 
Jn 


(2.24) 
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Assuming piecewise linear data, the fluctuation for the cell bounded by x* and Xj+i 
is evaluated as, 

9i,i + 1 = -A / u x di = -li i+ 1 = — AjF (2.25) 

./x, 

The elemental update, the LHS of Eqn. 2.24, is formed as, 

Snu t = £i ti+ 1 ^ j + «i+i t ) = 0i,i+i (2.26) 

Partitioning the fluctuation into halves and distributing equally to the nodes yields 
the basic non-upwind elemental update formula, 




^i.i + 1 4*1,1 + 1 

—Ui+lt - 


(2.27) 


Assembling all the elemental contributions to the nodal updates, it is clear each 
interior node will receive fluctuation signals from the elements adjacent to the left 
and right. The nodal update is formed as the sum of these fluctuation contributions, 


ti— l,i T Ii,i + 1 


^ it — S)Ui t — 


or, 

SiU it = + (2.29) 

A popular nomenclature convention for Eqns. 2.27 and 2.29 is to describe the elemen- 
tal distribution formula as, 

S iUit +- + COE, S i+lUi+u <- + COE (2.30) 

where COE indicates a sum of similar contributions from other elements joining at 
that node. 

Expanding the nodal update formula (Eqn. 2.29), 

o„. _ -ViF-AiF s ^ , 00 ^ 

Si u i i — 2 — (2.31) 

which is the identical central discretization as for DMFDSFV (Eqn. 2.10). 
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An upwind scheme can be constructed by introducing artificial dissipation in order 
to redistribute the fluctuation, 


e = sign(A)0 E 


(2.32) 


The upwind distribution formula becomes. 


SiUi t 


+ COE = LaiH - siell < A)) + COE 


Si+iu. 


i+lt 


e- 


4> E + <t> E 


+ COE = 


<t>i,i+ i(! + sign(A)) 


+ COE 


2 2 
Using the fluctuation definition (Eqn. 2.25) the nodal update is obtained as, 


(2.33) 


Si u it 


(A+|A|)VjU (A-lADAit* , p| |A| 2 _. 

— —Oit + —0;U 


(2.34) 


which is identical to the first-order upwind discretization for DMFDSFV (Eqn. 2.13). 

A second-order scheme is easily obtained by adding the exact same DMFDSFV 
correction, R 2 u (Eqn. 2.20), to the nodal update formula (Eqn. 2.34). 


2.2.2 Non-linear Advection 

Non-linear advection is obtained from Eqn. 2.6 by choosing the flux to be, 

F = y (2.35) 

Define the Jacobian of the flux, 

A = F u = u (2.36) 


so that, 


„ OF dF du „ 

F x "TT 7T F u U x 

ox du dx 


Equation 2.6 may be rearranged in non-conservation form, 


u t + F x = u t + Au x = 0 


(2.37) 
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DMFDSFV 

Following Roe[13], the analog to the numerical flux of Eqn. 2.15 becomes, 


fi 


i+i 


Ai+\A\i+\ A r — \A\ i+i 

2— «L + 2 UR 

F\_ + Fr l^li+7 , , 

(mr - ui) 


(2.38) 


2 2 

where A is the conservative linearization of the flux Jacobian, which in this case is, 

U\_ + Ur 


A+± 2 


(2.39) 


A first-order upwind scheme is obtained using piecewise-constant data, u\_ = u<, and 
ur = Ui + 1 . A second-order upwind scheme is constructed using the linear reconstruc- 
tion of Eqn. 2.16. The first-order residual may be written explicitly as, 


D X IP Mli-i 

R, = —6;F -l A,'/ V,;w 


2 


(2.40) 


Fluctuation splitting 

The elemental fluctuation is 

0E 


- / f x <m = - / 

Jn Jn 


Au r dD 


in Jq 

Assuming piecewise-linear data Eqn. 2.41 becomes, 

9e = -A i+ iAiU = -A iF 

An upwind scheme is created by introducing the artificial dissipation, 
4>'e = sign(i i+ i)^ E = -\A\ i+ iAiU 
The distribution formula remains, 

_ + COE 


(2.41) 


(2.42) 


(2.43) 


S%Ui t 


A/ +1 u-i + 1 t 


4>e + ( 1>'e 


+ COE 


(2.44) 
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The nodal update is, 


c _ 0i-l,i + §i- \,i , ( t > i,i + 1 
^ i u it o ' o 


V,F \A\ i_i 


V;U — 


A ,F |i| i+ i 


= SiF ~^V t u + 


2 

\A\ i+i 


+ 


A,u 


A ;u 


(2.45) 


This is the identical update formula as for DMFDSFV (Eqn. 2.40). 


Expansion shocks 


The discretization of Roe's scheme allows for non-physical expansion shocks that 
violate the entropy condition. Harten and Hyman [82] proposed a commonly used 
method for perturbing the wavespeeds such that entropy is satisfied and expansion 
shocks are prevented. The correction is applied to any wavespeed that can go to zero 
at a sonic point and takes the form, 


l^li+4 ^ 4 



where the perturbation scale is, 


if |A|j + i > e 
if |A| j+ i < e 


e 


= max 


0, 


( A *+i 


A,:), 



(2.46) 


(2.47) 


2.3 Scalar Advection-Diffusion 

The governing equation for scalar advection-diffusion problems in one-dimension is, 

u t + F x = {/ju x ) x (2.48) 


2.3.1 Heat Equation 

Modeling of the viscous RHS in Eqn. 2.48 begins with a consideration of the heat 
equation, 


u t = {im x ) x 


(2.49) 
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In the finite volume framework one approach to discretizing the viscous term is to 
construct a viscous flux, so that the nodal update becomes, 


SiU H = Q uu x ) i+ i - (im x ) t i 


(2.50) 


where, 


{T u x)i+i — Li + i 


(V«)i + (Vu) i+ i 


(2.51) 


with the gradients Vu defined by Eqn. 2.17. This approach leads to a five-point 
stencil. 

An alternative is to use a finite element discretization, which results in a three- 
point stencil. This approach is adopted both by Barth [22] and Anderson and Bon- 
haus[83] in a finite volume context and by Tomaich[50] in a fluctuation splitting 
context. 

A Galerkin finite element discretization, using mass lumping to the nodes, is 
constructed on the fluctuation splitting domain by integrating with the aid of the 
finite element linear shape function v (see Bickford §4.2.2[84] or Bathe §7.2[85]), 


S,u H 


/ Vi(im x ) t 

Jn 


dtt 


Integrating by parts, 


SiUi t 


= - / (vi) x (/m x ) dD 

J n 


(2.52) 


(2.53) 


The shape function is the linear tent function, and is equal to zero at Xi± i, eliminating 
the first RHS term of Eqn. 2.53. The remaining term is integrated over each element 
connecting at node i, 


‘d'Cii 


?/. 


v x u x fi dil 


(2.54) 


The dependent variable and shape function gradients are constant over the element, 
and taking the element-average viscosity coefficient the elemental contributions are, 

NjU 

2 


SiU it 

Si+l‘Ui+l t 


(-i,i + 1 

A jU 


l i,i + 1 


l*i+\ + COE 
Pi+- T COE 


(2.55) 
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The nodal update is written, 

SiUi t 


A, jU 
@i,i + 1 


A+§ 



1 _ 

2 


(2.56) 


2.3.2 Combined Advection and Diffusion 

The combined effects of advection and diffusion in the governing equation (Eqn. 2.48) 
are treated by discretizing the advection terms as discussed in section 2.2 and adding 
the discretization of the diffusion terms of section 2.3.1. Recall, however, that the 
upwind advection discretization includes artificial dissipation, which can mask the 
physical dissipation. 

Perhaps the best approach for solving discretized advection-diffusion problems, as 
suggested by Barth[21], is to include the maximum of either the physical diffusion 
term, as defined by Eqn. 2.51 or Eqn. 2.55, or the artificial dissipation, the second 
term of Eqn. 2.38 for DMFDSFV or A in Eqn. 2.43 for fluctuation splitting. 


2.4 Systems 

A hyperbolic conservation law for systems (Eqn. 2.2) is written in one dimension as, 

Uf + F iT =0 (2.57) 

A decomposition of the flux function is sought such that the system can be expressed 
as a decoupled set of advection-diffusion equations. 

2.4.1 Euler Equations 

The one-dimensional Euler equations[86] for perfect gases, suitable for simulating 
non-reacting, low-Knudson-number shock-tube flows, are written as a conservation 
law (Eqn. 2.57) with, 


U 


P 

pu 

pE 


> 


(2.58) 
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(2.59) 


pu 

F = <( pu 2 + P 
puH 

The Euler equations have a form similar to the non-linear advection problem. 

The total energy and enthalpy are obtained from the internal energy and enthalpy, 


u 


E = e+y 

The energy and enthalpy are related as, 


u 


h = e + — 
P 


H = h + - 


P 


The perfect gas equation of state is, 


p = pe(ry - 1) 


(2.60) 


DMFDSFV 

The numerical flux remains as in Eqn. 2.38, 

F|_ + Fr 


f, 


l+i 


(U R -U L 


(2.61) 


Roe[13, 80] constructs the conservative linearization for the |A| i+ i matrix by intro- 
ducing the parameter vector, 

I 1 ! 

Z = y/p ^ u (2.62) 

l H . 

The i + i state is taken to be a linear average of the parameter vector, 


z »+i = 


Z, — Z F 


Taking the velocity and total enthalpy from the parameter vector, 

~ 

«■=-=- H = 

Z\ 7i\ 


(2.63) 
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and defining the Roe-density, 

P = VplPr 

the Jacobian matrix is formed as, 

|A| = X|A|X -1 


The eigenvalues are, 

A = diag(u, u + a , u — a) 

The right eigenvectors are, 



f \ 

1 


( \ 

1 


1 

x (1) = < 

u 

> x (2) = < 

u + a 

> X (3) = < 

u — a 


u 2 

l 2 J 


H + ua 

V / 


H — ua 


The product X ^Ur — U L ) results in the characteristic variables, 


X _ 1 (U R — U L ) =x- ] dv = ~ 


2 a 2 dp — 2 dP 
dP + pa du > 
dP — pa du 


The sound speed is, 

a 2 = — = 7(7 - l)e = (7 - 1)^ = (7 -!)(#- \) 
P 2 

Also note the grouping pd,u can be constructed as, 


pdu = Z 1 dZ‘> - Z 2 dZ l 


(2.64) 

(2.65) 

( 2 . 66 ) 

(2.67) 

( 2 . 68 ) 

(2.69) 

(2.70) 


As for the scalar case, first-order spatial accuracy is obtained by taking the right 
state to be i + 1 and the left state at i. Higher-order accuracy is obtained using 
gradient reconstruction (Eqn. 2.16) applied either to each of the conserved variables 
(Eqn.2.58) or each of the primitive variables, which are, 


P 


V = < u 


P 


(2.71) 


The nodal update is still formed as in Eqn. 2.8. The residual remains as expressed 
in Eqn. 2.40, but for systems rather than scalar quantities. 
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Fluctuation splitting 

The Euler flux (Eqn. 2.59) can be written in terms of the parameter vector, 

F = } 2^ Zl Z 3 + ^-Zi (2.72) 

Further, the derivative of the flux is, 

Z 2 Z\ 0 

2±i Z 2 Zl dZ (2.73) 

0 Z 3 Z ‘2 

By assuming a linear variation of the parameter vector on each element, the 
fluctuation is obtained from Eqn. 2.41 as, 

(b e = — / Fj, dll = — I F zZ x dD — — F^AjZ (2.74) 

Jn Jn 

Deconinck[36] et al. show, 

F z AiZ = AA ( U = A.F (2.75) 

when the Roe-averaged forms (Eqns. 2.63 and 2.64) are used to obtain A. 

An upwind scheme is constructed by adding the artificial dissipation, 

4k = -|A|. + iA,:U (2.76) 

where |A| is defined in Eqn. 2.65. Employing the same distribution formula as for 
the scalar advection (Eqn. 2.44) leads to an update formula analogous to Eqn. 2.45, 
showing the equivalence between DMFDSFV and fluctuation splitting for the one- 
dimensional Euler equations. 

Before ending the fluctuation splitting discussion, it is desired to frame the arti- 
ficial dissipation in the form, 



</>'e = sign(A i+ i)<^> E 


(2.77) 
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The difficulty lies in defining the matrix sign(A). One approach combines Eqns. 2.65, 
2.74, 2.75, 2.76, and 2.77 to form, 

sign(A)A =|A| = X^IX” 1 

sign(A) = X|A|X _1 A _1 = XlAlA” 1 ^- 1 (2.78) 

Sidilkover[l] offers an alternative to brute force matrix multiplications for evalu- 
ating Eqn. 2.78. Introduce the auxiliary variables, W, defined by the transformation, 


dXJ = U wdW 


where, 


dYV 



with the first Riemann variable defined as, 


ds = dp — < -~ 


The Jacobian of the transformation is. 


Uw = 


1 0 

u 1 

,.2 


o" U 


1 i_ ft 

y— 1 “ r 2ft 2 


urv 1 = 


and its inverse is, 

(7-U* -(7-1)* 
—u 1 0 

(7-1 Yf -(7-1 )« 7-1 

The element fluctuation (Eqn. 2.74) can be reworked, 

<t> e = — AA,U = -UwU^AUwU^AiU = -U^AA,W = 
where 0 e is the fluctuation as computed for the auxiliary variables, 

= 


(2.79) 

(2.80) 


(2.81) 


(2.82) 


(2.83) 


(2.84) 


AAiW 


(2.85) 
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The flux Jacobian of the auxiliary variable formulation is obtained from the conserved 
flux Jacobian via the similarity transformation, 

A = U^AU w = Un/XAX-'U^ = XAX^ 1 (2.86) 

so the eigenvalue matrix, A (Eqn. 2.66), remains unaltered. The right eigenvectors 
are obtained from Eqns. 2.67 and 2.83, 


X = 


u^x 




The inverse is easily computed to be, 


\ 

1 


\ 

0 


/ 

0 

0 

> x (2) = < 

a 

> X^ = < 

—a 

0 

/ 


a 2 

\ > 


a 2 


X 


-i 


1 0 0 

n _L J_ 

U 2 a 2 a 2 

o __L J_ 

u 2 a 2 a 2 


The flux Jacobian is evaluated from Eqn. 2.86, 


A 


« 0 0 
0 u 1 

0 a 2 u 


(2.87) 


( 2 . 88 ) 


(2.89) 


which corresponds to the following non-conservative form of the Euler equations, 


s t + us x = 0 
pu t + upu x + P x =0 
P t + a 2 pu x + uP x = 0 


(2.90) 


Having developed an alternative method for obtaining the elemental fluctuation 
(Eqn. 2.84), the artificial dissipation can be addressed (Eqn. 2.77). 


b' E = sign (A) = U»-U M 1 sign(A)UiE^E = U wVi 


(2.91) 
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where, 

0e = U^ L sign(A)UvF<^ E = sign (A) (^ E (2.92) 

and with the aid of Eqns. 2.78 and 2.86, 

sign(A) =U^sign(A)U w = U^ 1 X|A|A- 1 X- 1 U„, = i’|A|A- 1 i’- 1 (2.93) 

Using the eigenvalue and eigenvector definitions (Eqns. 2.66, 2.87, and 2.88) sign(A) 
is evaluated to be, 


sign(A) = 


sign(d) 0 0 

0 | [sign('fi+a) + sign(d — d)] A [sign(d-fa) — sign (u — a)] 

0 | [sign(u+d) - sign(d — a)} \ [sign (ft +d) + sign(d — d)] 


(2.94) 


By considering two cases, for subsonic and supersonic conditions, Eqn. 2.94 takes on 
simple forms, 


\ i 

d| 

[ M sub i 

|u| 


where, 


and, 


sign (A) = 

M s “ p = sign(d)/ 


M s “ 6 = 


sign(d) 0 0 

0 0 \ 

a 

0 dO 


(2.95) 

(2.96) 

(2.97) 


2.4.2 Navier-Stokes Equations 

The Navier-Stokes equations[87, 88] for the flow of a perfect gas are written in one- 
dimensional conservation law form (Eqn. 2.57) with U defined in Eqn. 2.58 and the 
ffux defined as, 

F = F ?: — F" (2.98) 

where the inviscid flux, F z , is the same as the Euler flux (Eqn. 2.59). The viscous 
flux is, 

f 

0 


7~XX 

l U j 7~xx Qx 


(2.99) 
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Using Stokes’ hypothesis the stress is, 


Dx — g 


( 2 . 100 ) 


Fourier’s law for heat flow gives, 


Qx k,T x 


( 2 . 101 ) 


The thermal conductivity is related to the viscosity through the Prandtl number, 


( 2 . 102 ) 


where for air P r = 0.72[89]. The temperature is obtained from the perfect gas equation 
of state, 

T = A (2.103) 


The inviscid flux is discretized in the manner of section 2.4.1. The contributions 
from the viscous flux to the nodal update is obtained in a Galerkin sense using 
the system analog to Eqn. 2.54. No viscous contribution is made to the continuity 
equation. 

Using the linear variation of the parameter vector over an element, the velocity 
gradient is locally defined on an element [36, 34], 


where, 


£2 \ ___ z 2x z 2 , 

I 2 ^x 

ZlJ x Z X zf 


P = Z i 


T- {\z 2 - uAiZi) = A\u = iu (2.104) 
Zl zf p 


\fPi + \/Pi + 1 
2 


(2.105) 


is called the consistent density average. The viscous contribution to the momentum 
equation can now be expressed, 

r 4 4 _p 

/ -v x p,u x dD = -v x fi^A t u (2.106) 

J e 3 6 p 

The first term of the viscous energy flux is evaluated in a similar manner, 

/' 4 4 _ p 

/ - VxiJ,uu x di1 = -v x p,u—&iU (2.107) 

J e 3 3 p 
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The second term requires some manipulation. Begin by defining the temperature 
gradient over an element, 


(T x ) e = 


P 


X Px 

(M\ p ~p J 7 SR V p 


a 2 f AiP A ip 


P 


where, 


P = 


~2 ~ 

a p 

7 


The heat flow contribution to the viscous flux is then obtained, 


L 


v x kT x dil = v x k — 


a 2 ( AiP A ^ 


pc p a 2 ( AiP Aip 


P 


Pr 7 S J£ V P 


(2.108) 


(2.109) 


( 2 . 110 ) 


The elemental contributions from the viscous terms is similar to Eqn. 2.55, 


SiV k 




^i,i+ 1 


t- - 
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|u4Aj , u + ~p~ 
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V Si ( AiP _ 
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ti,i + 1 


0 

Ha iU 

3 D ' 


4 ~ P A 


Ct. 


+ COE 


+ COE (2.111) 


As discussed for the scalar advection-diffusion equations, when solving the Navier- 
Stokes equations the maximum of the viscous contribution to the nodal update and 
the artificial dissipation from the inviscid flux discretization should be utilized. When 
the physical viscous terms are large enough, no artificial dissipation is needed. 


2.5 Finite Volume State of the Art 

The one-dimensional analysis of the Navier-Stokes equations represents the state of art 
for upwind flux difference split finite volume schemes. Extensions of the unstructured 
finite volume method to multiple spatial dimensions relies upon solving a locally 
one-dimensional approximate Riemann problem across cell faces, an approach that 
looses some of the coupling present in a system of equations. Locally one-dimensional 
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solution techniques also introduce preferential grid-aligned wave directions that may 
not correspond with physical wave-propagation directions. 

The fluctuation splitting framework approaches the governing equations from a 
different perspective than finite volume, but is seen to result in identical discretiza- 
tions in one dimension for the DMFDSFV scheme. However, the following chapters 
show how fluctuation splitting generalizes to multiple dimensions in a more compact 
and coupled manner than locally one-dimensional finite volume schemes. 



Chapter 3 

Two-Dimensional Scalar Analysis 


Having shown the equivalence of the fluctuation splitting and DMFDSFV schemes 
for one-dimensional domains, the extensions to two spatial dimensions are considered. 
The present chapter analyzes the case for a single governing conservation law, 


Ut + V-F = V • (yuW) 


(3.1) 


to which steady-state solutions are sought. Systems, e.g. the Navier-Stokes equations, 
are deferred to a subsequent chapter. 

DMFDSFV is extended in an upwind, edge-based formulation for general unstruc- 
tured meshes with a multi-dimensional reconstruction. The crucial piece of the solver 
remains a locally one-dimensional approximate Riemann evaluator, as developed in 
chapter 2. This locally one-dimensional treatment of the fluxes results in increased 
production of artificial dissipation, particularly when discontinuities are not aligned 
with the mesh [40]. 

The extension of the fluctuation splitting scheme to multiple dimensions takes 
on the flavor of a node-based upwind residual-distribution algorithm, resulting in a 
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greater flexibility to propagate multi-dimensional wave phenomenon without dissipa- 
tion. Fluctuation splitting has a more-compact stencil than DMFDSFV for second- 
order accuracy and exhibits “zero cross-diffusion” 1 in a grid-aligned condition. Fluc- 
tuation splitting is seen to grid-resolve advection-diffusion problems on coarser meshes 
than DMFDSFV. 

This chapter begins by defining the different elemental domain structures for 
DMFDSFV and fluctuation splitting. Then formulations are developed for the two 
schemes and applied to two-dimensional linear and non-linear advection model prob- 
lems. Observations are made about the effects in both schemes of grid orientation 
on the production of artificial dissipation. Discretizations for diffusion follow, along 
with results for the heat equation. The chapter concludes with an advection-diffusion 
test problem, revealing the improved accuracy and decreased solution times using 
fluctuation splitting vis a vis DMFDSFV. 

3.1 Domain 

The DMFDSFV scheme is implemented for an edge-based data structure. The domain 
is discretized on an unstructured mesh of arbitrary connectivity. Control volumes are 
then constructed about each node. One common method for defining the control 
volumes is to use the median-dual mesh, shown as the dashed lines in Figure 3.1. 
For a triangulated domain, the generalized median-dual volume about a node equals 
one-third the sum of the areas of each triangle connected at that node, 

S, = l E S t (3.2) 

° vt | ie T 

The fluxes into and out of the control volumes are efficiently computed as a sum 
of contributions distributed to the nodes from a loop over edges. For each edge, 
the fluxes through the control faces to the right and to the left of the edge, Fig- 
ure 3.2, are computed, with the convention that a positive flux is out of the control 

1 “Zero cross-diffusion” refers to the practice of adding artificial dissipation terms in the stream- 
wise direction only, as opposed to adding artificial dissipation in both the streamwise and cross- 
stream directions. 
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(a) Interior node. 


• Node 

— Unstructured grid 

Median-dual control volume 

— ►- Outward normal 
x Quadrature point 

0 Current node 

1 Distance-one neighbor node 

2 Distance-two neighbor node 

interior 



Figure 3.1: Finite volume computational domain for edge-based implementation. 


volume surrounding the initiation node and into the control volume surrounding the 
termination node. The fluxes are evaluated at the quadrature points, and data is re- 
constructed from the nodes to the quadrature points, denoted by x , along the vectors 
f in Figure 3.2. 

For the special case of a boundary edge, shown in Figure 3.2(b), two fluxes are 
computed to the right-hand side of the edge, one for each associated node. 

The DMFDSFV scheme is referred to as a locally one-dimensional scheme because 
the fundamental Riemann problem is approximately solved normal to a face, with the 
solution going to either or both of two nodes only, connected by the physical mesh 
edge. The mesh edge and control volume face can have an arbitrary orientation within 
multi-dimensional space. The reconstruction step in general can be multi-dimension- 
al. 

In contrast to the edge being the fundamental computational element for DM- 
FDSFV, fluctuation splitting is formulated on a multi-dimensional simplex element. 
In two dimensions the simplex element is the triangle, while in three-dimensions the 
simplex is a tetrahedron. A pictorial of the domain nomenclature for fluctuation 
splitting is presented in Figure 3.3. Local curvilinear coordinates, are defined 
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Figure 3.2: Flux quadrature for edge-based finite volume scheme. Solid bold line is 
physical mesh, dashed lines are control-volume faces. 



Figure 3.3: Elemental triangular domain for fluctuation splitting. 
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on each triangle, parallel to two of the sides. A good choice of sides for the curvilinear 
coordinates to minimize computer round-off error may be the two most orthogonal 
sides. The fluctuation computed on a triangle may now be sent to one or more of 
the three vertices (versus two for DMFDSFV), allowing for a truly multi-dimensional 
distribution scheme. 

By definition, all elements are interior to the domain, so no special domain dis- 
cretization is needed for the boundaries with fluctuation splitting. However, numerical 
boundary conditions will still need to be applied at the boundary nodes to account 
for contributions from ‘ghost’ elements outside the computational domain. 


3.2 Advection 

Pure advection is obtained from Eqn. 3.1 when /< = 0, 

U t + V • F = 0 (3.3) 

This section extends the DMFDSFV procedure in a straight-forward manner from 
the one-dimensional analysis in chapter 2. Then the fluctuation splitting method is 
applied, in a formulation that now T differs significantly from DMFDSFV. A tempo- 
ral pseudo-time marching solution procedure follows, including a positivity analysis 
yielding timestep restrictions. A statement on the boundary conditions concludes the 
analytic formulations, leading to results for both linear and non-linear test cases. 

3.2.1 Formulations 
DMFDSFV 

The traditional, locally one-dimensional, approximate Riemann solver finite volume 
scheme[22] begins by integrating Eqn. 3.3 over the control volumes and applying the 
divergence theorem, 



(3.4) 
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Using mass lumping to the nodes, similar to an explicit finite element treat- 
ment [85], the temporal evolution is evaluated on a time-invariant mesh as, 

j[ U, da = s/IL -+ |r - Uf) (3.5) 

The discretization of the convective flux, F, is performed using Barth's imple- 
mentation^] of the upwind, locally one-dimensional, approximate Riemann solver 
of Roe[13] by constructing the numerical fluxes as a combination of the flux function 
and artificial dissipation, 


J F -hdT ~ y: / /ace AT = -R, (3.6) 

^ faces 

The numerical flux at the face is analogous to the one dimensional form (see Eqns. 2.12 
and 2.38), 

f face — Ty (^F in T F ou f''j ■ fl T (3. /) 

where the artificial dissipation provides the upwinding (see Eqns. 2.14 and 2.38), 

T = l -\X-h\(U 0 ut-U m ) (3.8) 

Out and in refer to states on the outside and inside of ii at the face. The flux 


Jacobians are defined, 


dF x 8F y 

A* = , ,4 y = 


dU ’ ‘ dU v J 

where A = AO + A y j, F = F x i + F y j, and the tilde indicates the conservative 
linearizations at the cell face [13]. 

Piecewise linear reconstruction from the nodal unknowns to the cell faces as in 


Eqn. 2.16, repeated here, 

Uface = U t + ipNU ■ r (3.10) 


provides second-order spatial accuracy in smoothly-varying regions of the solution. 
Median-dual gradients of the dependent variable, VC/, are obtained from the un- 
weighted least squares procedure outlined by Barth. Following Bruner and Wal- 
ters[90], the limiter is supplied an argument equal to half the argument Barth uses, 


namely, 



(3.11) 
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where U mtn l max is the minimum (resp. maximum) of Ui and all distance-one neighbors. 
The most restrictive limiting from choosing the minimum or maximum is used. 

In casting the limiter argument in this form, Bruner equates the Barth limiter 
with Superbee, for a limiter argument less than or equal to one. For the full domain 
of the argument, the non-symmetric Barth limiter takes the form, 

{ 0 | < 0 

2£ if 0 < £ < i (3.12) 

Q ~ 2 

for the limiter cast as Eqn. 3.11. Many limiter functions exist, and several of the 
more popular versions are detailed in appendix A. 

The DMFDSFV flux evaluation needs to be performed three times for each inte- 
rior triangle, once for each edge of the triangle. For linear advection at a uniform 
advection speed Eqn. 3.8 requires 4 multiplication/divisions and 2 additions/sub- 
tractions. Equation 3.7 requires 4 each multiplications and additions. Ignoring the 
work required to compute the nodal gradients and limiter, which varies based upon 
the mesh connectivity, the reconstruction, Eqn. 3.10, requires 6 multiplications and 
4 additions per face. Equation 3.6 adds one more multiplication per face, bringing 
the total operation count for DMFDSFV per triangle to 45 multiplications and 30 
additions. 


Fluctuation splitting 


The Narrow Non-Linear (NNL) fluctuation splitting scheme is presented as a slight 
re-interpretation of the work of Sidilkover and Roe[3]. The current interpretation 
is as a volume integral over triangular elements, without recourse to the divergence 
theorem. The discretized equations, however, are identical to Sidilkover’s. This form 
of fluctuation splitting employs a general limiter function for determining the residual 
distributions. 

Integrating Eqn. 3.3 over an element, where D is now the area of the triangular 
element, 


/ U t dtt =- / 


v ■ f dii 


(3.13) 
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For linear variation of the dependent variable over the element, the temporal evolution 
is (see also Eqn. 2.26), 

/ U t dtt = SrU t = ^-(U h + U 2t + U 3t ) (3.14) 

Jn 3 

where Ui, U 2 , and U 3 correspond to the three nodes defining element O. 

Using the local curvilinear coordinates (£, 77 ), defined in Figure 3.3, the divergence 
of the convective flux can be transformed, 

V . F = F; + FI = A ( p . i) + |-(F . j) (3.15) 


A ] 

dx l _ 

A [ 

.dy) 


tr, Vx 

'Ii 


d_ 

d_ 

. dr] . 


V - F = £ X (F ■ ?) f + £y{F ■ j) f + >lx(E ■ *)„ + Vy(F ■ j)y 
Introducing the inverse Jacobian of the coordinate transformation, 

t-i _ _ 25't 

■J — x ^lJrj 


and the invariants of the transformation, 

Vn _ hi A 


& = 
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-1 
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-1 
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Vi _ «3-v 


f X V n Vj 
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-1 
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Eqn. 3.17 becomes, 


Vy = 


V-F =j^(nr4-n 3 -4) =^(nrlt/ c -n 3 --4f/. 


J 


-1 


J- 1 


(3.16) 

(3.17) 

(3.18) 


(3.19) 


(3.20) 


If F is linear or quadratic in U, then for a linear variation of U over the element, 

I V -Fdtt = aA/rU + [3\U (3.21) 

Jn, 

where the difference operators are A (U = U- 2 — U\ and A V U = U 3 — U 2 and the 
advection speeds are, 


t m x 
a = —ni-A, 


o h ~ T 

P = ~~ n 3 -A 


(3.22) 
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A is now the conservative linearization over the entire triangular element [36]. 

The advective fluctuation can be defined, 

(j) =- f V -Fdtt (3.23) 

Jo. 

The fluctuation can be split. 

<f> = <f>* + <F (3.24) 

where, 

H = —aA^U, <F = -f}\U (3.25) 

Following Sidilkover[l] the scheme is extended to second-order spatial accuracy 
by repartitioning the fluctuation through the use of a symmetric limiter function, 

^ + o"r { Q) =: / ( 1 - (3.26) 

= r - <fhq) = r (i - hq)) ( 3 . 27 ) 


with, 



(3.28) 


In practice, if an averaging function, M,i , , exists for the desired limiter, it is numerically 
advantageous to compute M^(^>' ) , — H) = d r >'ip( Q ), avoiding the need to evaluate Q 
explicitly. 

This critical step, allowing the redistribution of the fluctuation, is what principally 
distinguishes the multi-dimensional fluctuation splitting scheme of Sidilkover from a 
locally one-dimensional extension of Riemann solvers. There is no analog to this in 
the formulations of chapter 2. 

Upwinding is achieved through the introduction of the artificial dissipation terms, 


<p h = sign (a)(f)*\ 4>'" = Agn(p)(f)*" (3.29) 


Combining Eqn. 3.14 with a distribution scheme for Eqn. 3.23 and summing over 
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all elements, the contributions to nodal time derivatives can be written in the form, 

s\r/ lt <- htf* -<ff) + coE 

S 2 U 2t <r- + ^) + \{cb* v -f) + COE 

SsU 3t <- + f) + COE (3.30) 

or in a more compact form, 

SiU H <- i [i(3 - i)(^ + (-1)'/) + (-4 + 5 * - i 2 )(<T - (-l)V'")' 

+ COE i =1,2,3 (3.31) 

where COE stands for contributions from other elements containing these nodes. 

The distribution requires 4 addition/subtractions and 3 multiplication/divisions. 
Upwinding requires 2 multiplications for each term of Eqn. 3.29. Ignoring the cost of 
evaluating the limiter function, as was done when counting the DMFDSFV operations, 
one more multiplication and addition are performed in each of Eqns. 3.26 and 3.27 and 
for each term of Eqn. 3.25. Finally, a and [L Eqn. 3.22, each require 4 multiplications 
and one addition for a total operation count of only 19 multiplications, versus 45 for 
the DMFDSFV scheme. Only 10 additions per triangle are required by fluctuation 
splitting, versus 30 for DMFDSFV. 

3.2.2 Advective Timestep Restriction 

Both schemes are formulated as either Gauss-Seidel or Jacobi time-relaxation algo- 
rithms. 

The nodal updates for the discrete system can be formed as a sum of contributions 
from all nodes. 

U‘ +M = V c i u i = “Vi + E c i U > (3.32) 

3 3 & 

For positivity[56] each of the coefficients in Eqn. 3.32 must be non-negative. 

In the finite volume context the nodal update (Eqn. 3.32) can be rearranged into 
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the form of Eqn. 3.5, 

§- t a t +i ‘ - u>) = ^(c - m + §- t E °i u i < 3 - 33 > 

For the upwind, edge-based algorithm considered here, each -A-Cj will be related to 
a positive-definite coefficient equal to zero for outflowing faces and related to the 
wavespeed for inflowing faces, yielding the restriction At > 0 on the timestep. The 
remaining term can be expressed, 

^(q- 1) =-^c fc (3.34) 

k about i 

where the c coefficients are also positive-definite, either zero for inflowing faces or re- 
lated to the wavespeed for outflowing faces. Rearranging and imposing the positivity 
constraint, c$ > 0, yields the timestep restriction, 

1 - — Yl Ck = C i A 0 (3-35) 

1 k about i 

At < (3.36) 

£—sk about i ^k 

For fluctuation splitting, the nodal updates are assembled from Eqn. 3.31 as, 

§(ur M -u>) = £>((/,- 1/,-) (3.37) 

In this case the Cj coefficients are formed as contributions from the fluctuations in 
the triangles to both the left and the right of mesh edge Vj. The positivity restriction 
on At is found to have a similar form as for finite volume (Eqn. 3.36), 

At < — ^ — (3.38) 

Local time-stepping based on positivity is shown to yield stable, yet non-conver- 
ging, solutions in some second-order cases (see section 3.4). Robust convergence is 
obtained by using the first-order c’s in Eqns. 3.36 and 3.38, even for second-order- 
accurate spatial discretizations. This is similar to the common practice of using a 
first-order Jacobian discretization in a time-implicit scheme. 



42 


CHAPTER 3. TWO-DIMENSIONAL SCALAR ANALYSIS 


An implicit scheme can be constructed for fluctuation splitting by linearizing the 
flux in time. 


pi+At = pt + (Jjt+At _ jji) ( 3 . 39 ) 

The spatial discretization remains the same and the coefficients of the LHS of Eqn. 3.37 
are changed from the scalar ^ to the operator, 

— + aA^ + /3A,] (3.40) 

A full implicit scheme would require inverting the sparse matrix corresponding to 
the operator of Eqn. 3.40. Since both and A v operators use a nearest-neighbor 
compact stencil, creative re-ordering of the nodes could pack the implicit matrix into 
a banded diagonal of average bandwidth 6. Another strategy is to use a variant of a 
red/black coloring, here grouping the nodes into three sets because of the triangular 
connectivity. The three sets of nodes would then be integrated with a line Gauss- 
Seidel strategy, where the matrix to invert has been reduced to a simple diagonal form 
for each set of nodes. The present dissertation does not implement a fully-implicit 
temporal integration. 

3.2.3 Boundary Conditions 

Explicit Dirichlet inflow boundary conditions are employed. Advective outflow bound- 
aries are treated for free convection through the boundary nodes, allowing boundary 
nodes to be handled in the same manner as interior nodes. 

3.2.4 Results 

Linear 

The linear advection equation is obtained from Eqn. 3.3 for a flux function F = XU, 
yielding, 


U t + V-(XU) =0 


(3.41) 
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Figure 3.4: First-order fluctuation splitting, uniform advection. 


A divergence-free advection velocity is considered, such that V-A = 0. Equation 3.41 
can then be written, 

U t + X-VU =0 (3.42) 

The appropriate averaged flux Jacobian for linear advection is simply A = A, 
where A is evaluated at the quadrature points for DMFDSFV and at the element 
centroid for fluctuation splitting. 

Uniform advection 

Uniform advection of the Heavyside function at —45 degrees, A = (1, —1), on a cut- 
cartesian mesh is shown for first-order fluctuation splitting, second-order fluctuation 
splitting, and second-order DMFDSFV in Figures 3. 4-3. 6, respectively. The mesh 
is shown as the dashed background, and equally-spaced contours vary on [0,1], the 
minimum and maximum solution values. The spread of the contour lines with spatial 
evolution is indicative of the amount of dissipation introduced into the solution by 
the discretization of the convective terms. 
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Figure 3.5: Second-order fluctuation splitting, uniform advection. 


Second-order fluctuation splitting is seen to be greatly superior to first-order, as 
expected, reproducing the exact solution in this case with no introduced dissipation. 
Also, fluctuation splitting represents a significant reduction in numerical diffusion ver- 
sus the corresponding DMFDSFV scheme, with both results employing the Minmod 
limiter. 

However, the “zero cross-diffusion” results of Figure 3.5 with fluctuation splitting 
are misleading. In Figure 3.7 the advection velocity has been rotated counter clock- 
wise by 90 degrees on the same grid. Clearly, the artificial dissipation introduced by 
the fluctuation splitting scheme has been increased. 

The corresponding DMFDSFV solution is shown in Figure 3.8. While the change 
in contour spreading for the DMFDSFV scheme between Figures 3.6 and 3.8 is less 
dramatic than the change in spreading for the fluctuation splitting scheme in Fig- 
ures 3.5 and 3.7, the fluctuation splitting results still exhibit less diffusion than the 
DMFDSFV results, comparing Figures 3.7 and 3.8. 

Employing the compressive Superbee limiter with the fluctuation splitting scheme 
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Figure 3.6: Second-order DMFDSFV, uniform advection. 


yields the results of Figure 3.9. In this case the discontinuity is confined to a 2-3 cell 
stencil, and does not grow in space. Applying the Superbee limiter to DMFDSFV 
cannot eliminate all artificial dissipation for this case, as is possible with fluctuation 
splitting. The DMFDSFV results (not shown) corresponding to Figure 3.9 spread the 
discontinuity over four cells by the outflow boundary, with a continually broadening 
trend. 

However, while it is possible to use the Superbee limiter with fluctuation splitting 
for this case, compressive limiters can be unstable on different grid orientations. For 
example, no degree of compression is stable for the case of Figure 3.5. This potential 
for instability is related to global positivity, as discussed by Sidilkover and Roe[3]. 

The effect of using a general unstructured grid is investigated in Figures 3.10 
and 3.11. The unstructured grid in this case was generated using Vgrid[91, 92]. 
The fluctuation splitting solution exhibits less dissipation, but is not as smooth as 
the DMFDSFV solution. While the fluctuation splitting scheme preserves contact 
discontinuities over larger spatial ranges than the DMFDSFV scheme, fluctuation 
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* = (1.1) 

Contour 

spacing: 

0.05 - 0.95 
0.1 increments 


Figure 3.7: Second-order fluctuation splitting, uniform advection. 



* = (i.i) 

Contour 

spacing: 

0.05 - 0.95 
0.1 increments 


Figure 3.8: Second-order DMFDSFV, uniform advection. 
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Figure 3.9: Second-order fluctuation splitting with compressive limiter. 


splitting does not appear to degenerate gracefully with regard to extreme coarsen- 
ing of the unstructured mesh for this test case. This behavior could have negative 
implications for applications employing multigrid convergence acceleration. 

Circular advection 

Circular advection is achieved by setting A = ( y , —x). A decaying sine- wave input 
profile is used, 

U (. x , 0) = (e x sin irx) 2 

Results for the two schemes, using the Minmod limiter, are presented on the worse- 
case cut-cartesian mesh in Figures 3.12 and 3.13. Again, the fluctuation splitting 
results are considerably less diffusive than the DMFDSFV solution. 

The circular-advection problem is also applied on an unstructured mesh. The 
input profile for this case consists of both a top-hat function and a decaying sine 
wave, allowing comparisons between the schemes for both sharp discontinuities and 
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* = (1.1) 

Contour 

spacing: 

0.05 - 0.95 
0.1 increments 


Figure 3.10: Fluctuation splitting on unstructured mesh. 



* = (i.i) 

Contour 

spacing: 

0.05 - 0.95 
0.1 increments 


Figure 3.11: DMFDSFV on unstructured mesh 
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X = ( y, -x ) 

Contour 
spacing: 
0.025-0.375 
0.05 increments 


Figure 3.12: Fluctuation splitting, circular advection. 



x = ( y, -x ) 


Contour 
spacing: 
0.025-0.375 
0.05 increments 


Figure 3.13: DMFDSFV, circular advection. 
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smooth gradients. The input profile is, 


U(x, 0) 


( (e 2x sin(27r.r)) 2 


< 


0 

0.4 

0 


-0.5 < x <0 
—0.6 < x < —0.5 
—0.8 < x < -0.6 
— 1 < x < —0.8 


Results for this case are displayed in Figure 3.14 for fluctuation splitting and Fig- 
ure 3.15 for DMFDSFV, both using the Minmod limiter. Fluctuation splitting per- 
forms significantly better at preserving the top-hat distribution. Fluctuation splitting 
also does a better job of maintaining the minimum and maximum values of the sine 
distribution, though both schemes do well on the smooth gradient portion of the sine 
wave. 


Non-linear 

The non-linear advection equation is obtained from Eqn. 3.3 for the flux function 
F = (—. U). In non-conservative form the equation is written, 

Uf, + UU X + Uy = 0 

A coalescing shock problem is considered, with an anti-symmetric input profile, 

U(-l,y) = U(0,y) = 0 
U(x, 0) = — 2x — 1 on x = (—1, 0) 

The exact solution to this problem contains symmetric expansion fans on the sides and 
a compression fan at the inflow that coalesces into a vertical shock at ( x , y) = (— ^). 

The first mesh is cut-cartesian containing 26 x 26 nodes. The fluctuation splitting 
and DMFDSFV solutions, both using the Minmod limiter, are presented in Fig- 
ures 3.16 and 3.17. Both algorithms exhibit the same grid dependence on the amount 
of artificial dissipation as seen before, with the left-half solutions having more diffu- 
sion than the right halves, due to the grid orientation. Both methods perform the 
same in the compression-fan region, coalescing into a shock to within the accuracy of 
the input-profile discretization. 
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X = ( y, -x ) 

Contour 
spacing: 
0.025-0.375 
0.05 increments 


Figure 3.14: Fluctuation splitting on unstructured mesh, circular advection. 


x = ( y, -x ) 

Contour 
spacing: 
0.025-0.375 
0.05 increments 



Figure 3.15: DMFDSFV on unstructured mesh, circular advection 




Figure 3.16: Fluctuation splitting, Burgers equation. 



Figure 3.17: DMFDSFV, Burgers equation. 
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The shock is more sharply defined by fluctuation splitting than by DMFDSFV. 
Figure 3.16 has the correct shock speed, with nearly the entire gradient captured in 
one cell thickness. In contrast, Figure 3.17 shows a slightly incorrect shock speed when 
using DMFDSFV, as the shock progresses to the left beyond the coalescence point, 
even though the discretization is conservative. The incorrect shock speed results from 
a non-symmetric distribution of the dependent variable to the left and right of the 
shock, caused by the excessive artificial diffusion generated on the grid-misaligned 
(left-hand) side. 

Contours of the absolute value of the error are presented in Figures 3.18 and 3.19. 
Errors from both computed solutions show a lack of symmetry, again reflecting the 
grid dependence of the artificial diffusion terms. The error levels from fluctuation 
splitting are less than from DMFDSFV. The shock curvature in the DMFDSFV 
solution at the coalescing point is clearly visible in Figure 3.19, resulting in significant 
downstream errors in the shock location as compared with the fluctuation splitting 
errors. 

This problem is repeated on a 25 x 25 mesh with symmetric diagonal cuts, favor- 
ably aligned with the advection directions. The fluctuation splitting and DMFDSFV 
solutions, Figures 3.20 and 3.21, are in good agreement. Plots of the absolute error 
contours, Figures 3.22 and 3.23, show fluctuation splitting to be a little more accurate 
than DMFDSFV for this case. 

The final mesh for this case is a truly unstructured triangulation containing 
847 nodes and 1617 cells. The nodes are clustered to the outflow boundary, with 
an arbitrary bias towards the left-hand side, introduced merely as an additional chal- 
lenge for the schemes. The fluctuation splitting solution is presented in Figure 3.24, 
showing very accurate and crisp shock resolution and good symmetry in the solution 
contours despite the mesh-clustering bias. In contrast, the DMFDSFV solution in 
Figure 3.25 has a more-diffuse shock and again an incorrect shock speed, with the 
outflow shock offset to the left of x = — The DMFDSFV solution is also somewhat 
less symmetric than the fluctuation splitting solution. 
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Figure 3.20: 
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Fluctuation splitting, Burgers equation, symmetric mesh. 
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-1 - 1 

0.1 increments 


U = 0 


Figure 3.21: DMFDSFV, Burgers equation, symmetric mesh. 
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Figure 3.24: Fluctuation splitting, Burgers equation, unstructured mesh. 



Figure 3.25: DMFDSFV, Burgers equation, unstructured mesh 
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3.3 Diffusion 

The model diffusion equation, the well known heat equation, is obtained from Eqn. 3.1 
when F = 0, 

u t =V-(/iW) (3.43) 

Lacking cross-derivative terms, the heat equation suffers somewhat as a model for the 
Navier-Stokes equations. 

Two methods for evaluating the diffusion term are detailed. The finite element 
form has a more compact stencil and will be shown to be more accurate. Either dif- 
fusion discretization can be used with either advection discretization, and in practice 
the finite element diffusion treatment is preferred with both the fluctuation splitting 
and DMFDSFV advection schemes. 

3.3.1 Formulations 

Finite volume 

In finite volume form, with the aid of the divergence theorem, Eqn. 3.43 becomes, 

/ U t dil = l ijNU ■ h flY (3.44) 

./o ,/r 

Discretizing the RHS of Eqn. 3.44, in a manner similar to Eqns. 3.6 and 3.7 (see also 
Eqns. 2.50 and 2.51) yields, 

f l^U ■ h dT A (? Uin + ^ Uoul ) ' " Ar (3-45) 

^ faces 

where in refers to the state inside the control volume and out refers to the state 
on the other side of the face from the control volume. The diffusion coefficient is 
averaged over the length of the face. The gradients from Eqn. 3.10 are not limited 
before averaging at the control-volume faces in Eqn. 3.45, as suggested by Anderson 
and Bonhaus[83]. 
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Finite element 


A finite element treatment, similar to Tomaich[50], is employed. Begin by integrating 
Eqn. 3.43 over an element, where if is now the area of the triangular element, 

/ U t dil = I V-(/iW)tft2 (3.46) 

Jn Jn 

The diffusive fluctuation is defined, 

<p v = / V • (fjNU) (in (3.47) 

Assuming piecewise-linear data and an element-averaged diffusion coefficient leads 
to a diffusive fluctuation of zero for the triangular element. Introducing the linear 
nodal shape functions such that i v i = 1; the elemental diffusive fluctuation 
can be expressed (ff = Vi = 0> where 

$ = / ViV ■ (ijVU) cin (3.48) 

Jn 

Integrating by parts (analogous to Eqn. 2.53), 


<f> v i = j vmVU -hclY - J nfU-VvidO (3.49) 

The boundary integral in Eqn. 3.49 will cancel for interior nodes on summing con- 
tributions from neighboring elements. For linear variation of U over the element the 
remaining volume integral can be evaluated analytically, 


w =D-^Ujt 


2S t 


« =7mt, u >- 


J n j 


Vu,- = 


—iihi 


j = i 


25 t 


3 = 1 


/."‘—fig 


Ujijhj-hi 


(3.50) 


(3.51) 


where //, is the element-average diffusion coefficient. 

The physical dissipation is then distributed to the nodal updates as, 


S t U H /— $ + COE 


(3.52) 
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3.3.2 Diffusive Timestep Restriction 


Continuing the positivity discussion from section 3.2.2, restrictions on the timestep 
are sought to maintain positive coefficients c* in Eqn. 3.32. 

Unfortunately, the finite element formulation for the diffusive terms (Eqn. 3.51) 
cannot be guaranteed to preserve local positivity on obtuse triangles (see Barth[22]). 
Considering only the contributions from the current node, the coefficient for the 
diffusion term can be written, 


Uj+At 



The resulting timestep restriction is, 


At < 


Si 

y- PL 

■ 'T 45t 


In a similar manner the timestep restriction from Eqn. 3.45 is, 


At < 


Sj 

3/i(Ar) 2 

4SV 


(3.53) 


(3.54) 


(3.55) 


3.3.3 Boundary Conditions 

For the diffusion terms a Neumann outflow boundary is applied with zero gradient 
normal to the boundary (VU-h = 0), achieved by setting the boundary integral in 
Eqn. 3.49 to zero. 


3.3.4 Heat Equation 

The heat equation (Eqn. 3.43) test problem, a steady-state boundary value problem 
on the unit square in the second quadrant, is taken from Tomaich[50]. The Dirichlet 
boundary values are, 


U(0,y) = sin (ny) 
U(x, 1) =— sin(7ra:) 


U(-l,y)=0 
U(x, 0) = 0 
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The analytical solution on x = [— 1) 0], y = [0, 1] is, 

U(x.y ) = — - — fsinh(7r(x+l)) sin(7ry) + sinh(7ry) sin(7r(x+l))l 
smli 7T 

Both diffusion discretizations, Eqns. 3.45 and 3.51, are compared on a 438-node 
unstructured mesh. Figures 3.26 and 3.27 plot the absolute value of the error in 
the converged solutions using Eqns. 3.45 and 3.51, respectively. A carpet plot of the 
solution, using the finite element formulation, is presented in Figure 3.28. 

The finite element treatment is clearly more accurate, and is used to discretize 
the diffusion terms for both DMFDSFV and fluctuation splitting in the following 
section. The average-gradient results in Figure 3.26 appear to exhibit a decoupling 
mode, similar to odd/even decoupling for structured meshes. 


3.4 Advection-DifFusion 


When considering combined advection and diffusion problems, rather than adding 
both physical and artificial dissipation terms, it is advantageous to only augment 
the physical dissipation with as much artificial dissipation as is necessary to ensure 
monotonicity of the solution. That is, if the physical dissipation term is larger than 
the computed artificial dissipation, no artificial dissipation is needed. 

This rationale can be implemented with fluctuation splitting, for example, by 
modifying the fluctuation distribution (Eqn. 3.30) to be, 


SiU lt 

S 2 U 2t 

s 3 u 3t 


^- + max^-^-,^j +COE 
^ + J* + max ( W 4>l )+ COE 


n 11 


+ max r <f>l ) + COE 


(3.56) 


Similarly for DMFDSFV, the artificial dissipation term in Eqn. 3.8 can be compared 
with the physical diffusion term. 
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Figure 3.26: Pure- diffusion problem error, diffusion terms from Eqn. 3.45. 
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Figure 3.27: Pure-diffusion problem error, diffusion terms from Eqn. 3.51. 
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0.75 


0.5 


Figure 3.28: Heat equation solution using finite element formulation. Contour incre- 
ment is 0.1. 

Smith & Hutton problem 

The final scalar test case is a linear advection-diffusion (Eqn. 3.1) problem due to 
Smith and Hutton [93]. The flux function is F = \U, with, 

A = (2f/(l - x 2 ), —2x(l - y 2 )) 

The streamlines for this problem, while not truly circular, are similar in orientation 
to the circular advection problem. The inflow profile is, 

U{x, 0) = 1 + t.anh(20x + 10) 

The diffusion coefficient is chosen to be a constant, y, = 10 -3 . The domain is the unit 
square in the second quadrant. No closed-form solution is known to the author for 
this problem. 

A sequence of five unstructured meshes is considered. The meshes have no pre- 
ferred clustering or stretching and have nominal node-spacings of 0.1, 0.05, 0.025, 



64 


CHAPTER 3. TWO-DIMENSIONAL SCALAR ANALYSIS 


Mesh 


A 

B 

C 

D 

E 


Nodes Fluctuation Splitting DMFDSFV 
(CPU seconds) 


134 <1 <1 

495 1 1 

1,928 5 8 

7,529 64 145 

28,915 760 1880 


Table 3.1: Grids and solution times for advection-diffusion problem. 


Fluctuation Splitting DMFDSFV 

\n 11^11-2 Mesh ||*|| 2 ||<f|| 2 

(art.) (phys.) (art.) (phys.) 


1274 

215 

A 

1918 

190 

597 

265 

B 

640 

176 

192 

161 

C 

144 

119 

54 

76 

D 

46 

66 

13 

36 

E 

18 

36 


Table 3.2: L 2 -norms (xlO 5 ) of artificial and physical viscosities for advection-diffusion 
problem. 


0.0125, and 0.00625, labeled as Meshes A-E, respectively. The number of nodes for 
each mesh, along with the solution times for both fluctuation splitting and DMFDSFV 
on a 195 MHz SGI R10000 CPU are listed in Table 3.1. 

L 2 -norms of the artificial and physical viscosities computed using both fluctuation 
splitting and DMFDSFV are presented for each mesh in Table 3.2. Notice that the 
norm of the artificial dissipation for both DMFDSFV and fluctuation splitting drops 
lower than the norm of the physical dissipation on Meshes D and E. However, the 
norm of the physical dissipation is smaller for DMFDSFV than fluctuation splitting on 
each mesh A-D. The physical viscosity is driven by the solution curvature, suggesting 
fluctuation splitting maintains the solution profile sharper than DMFDSFV" on the 
coarser meshes. A comparison of outflow profiles will soon verify this interpretation. 

Evidence of a grid-resolved fluctuation splitting solution is seen in Figures 3.29 
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Figure 3.29: Fluctuation splitting profiles on finest mesh, advection-diffusion problem. 
(FS = fluctuation splitting) 


and 3.30. The fluctuation splitting solution on Mesh E at the outflow boundary is 
presented along with the inflow profile and the corresponding pure-advection (/j = 0) 
fluctuation splitting solution in Figure 3.29. Because the pure-advection solution 
is seen to replicate the inflow profile, the fluctuation splitting artificial dissipation is 
insignificant on this mesh, and hence further grid refinement would not appreciably 
change the solution. Plotting only the fluctuation splitting results with respect to 
grid refinement, Figure 3.30 shows a convergence of the outflow profile by Mesh C for 
fluctuation splitting. 

The accuracy of fluctuation splitting and DMFDSFV are compared in Figure 3.31, 
where the outflow solutions from fluctuation splitting and DMFDSFV are plotted for 
Meshes C and E. Taking the grid- converged fluctuation splitting Mesh-E solution to 
be the true solution, it is clear that fluctuation splitting reaches the grid converged 
solution on a coarser mesh than DMFDSFV. 

Computational efficiencies of the two algorithms are compared in Figure 3.32, 
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Figure 3.30: Fluctuation splitting grid convergence, advection-diffusion problem. 



Figure 3.31: Fluctuation splitting and DMFDSFV for advection-diffusion problem. 
(FS = fluctuation splitting, FV = finite volume) 
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Figure 3.32: Convergence histories for advection-diffusion problem. (FS = fluctuation 
splitting, FV = finite volume) 


where the Z, 2 -norm of the residual is plotted versus CPU time for the fine-mesh 
fluctuation splitting and DMFDSFV solutions, along with the fluctuation splitting 
convergence history on Mesh D. The Mesh-E fluctuation splitting solution converges 
in 760 sec. The corresponding DMFDSFV solution takes 2.5 times longer than fluc- 
tuation splitting, due, in part, to the need to reconstruct gradient information at 
each node with DMFDSFV for second-order spatial accuracy. However, considering 
the solution time to reach a given accuracy, it is more reasonable to compare the 
fluctuation splitting solution time on Mesh D to the finest-mesh DMFDSFV solu- 
tion. The fluctuation splitting Mesh-D solution took only 64 sec, a factor of 29 times 
less than DMFDSFV on Mesh E, and still shows better accuracy than the fine-mesh 
DMFDSFV solution. 

An even greater speedup is seen with fluctuation splitting in conjunction with the 
Van Albada limiter, where now the Mesh-B solution over-plots the curve from the 
finest grid, shown in Figure 3.33. The corresponding DMFDSFV result using the 
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Figure 3.33: Advection-diffusion results using Van Albada limiter. (FS = fluctuation 
splitting, FV = finite volume) 


Van Albada limiter on Mesh B is included, and clearly falls short of the fluctuation 
splitting accuracy. The DMFDSFV case was repeated with the highly-compressive 
Superbee limiter with little improvement in accuracy. The solution time for fluctua- 
tion splitting on Mesh B is about one second, yielding a speedup factor of 2-3 orders 
of magnitude over DMFDSFV. 

The final set of results addresses convergence issues while pushing the positivity 
limits. Figure 3.34 compares two convergence histories for the second-order fluctu- 
ation splitting on Mesh B. The non-converging, though stable, convergence history 
is the result of using strict positivity arguments to set the timestep (Eqn. 3.38). 
The resulting solution is bounded and approximately correct but oscillatory. Limiter 
“ringing” is the cause of this behavior, with the higher-order discretization for the 
implicit matrix reducing the diagonal dominance, and hence stability, of the Gauss- 
Seidel iteration. 

Full convergence is achieved by using first-order positivity coefficients, which are 
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Figure 3.34: Convergence rates using first- and second-order positivity coefficients. 


not dependent on the limiters. The resulting local timesteps will not be as large as 
true second-order positivity would allow, but are more robust. 


3.5 Benefits of Fluctuation Splitting 

Fluctuation splitting and DMFDSFV schemes have been compared in detail as ap- 
plied to scalar advection, diffusion, and advection-diffusion problems. The fluctuation 
splitting scheme is seen to introduce less artificial dissipation while treating advection 
terms, allowing for more accurate resolution of weakly dissipative advection-diffusion 
problems. The ability to resolve solutions to these problems on coarser meshes makes 
the fluctuation splitting scheme the preferred choice over DMFDSFV. 

Linear advection test problems were utilized to investigate the dependence of 
artificial diffusion production on grid orientation. Both fluctuation splitting and 
DMFDSFV are shown to exhibit grid dependencies, but with fluctuation splitting 
producing less artificial dissipation on all grids considered. 
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A non-linear coalescing shock problem further explored grid dependencies as cases 
were constructed that result in incorrect shock speeds for DMFDSFV. Fluctuation 
splitting shows correct shock speeds for all grids and provides tighter shock capturing 
than DMFDSFV. 

An advection-diffusion problem with small physical dissipation (diffusion coeffi- 
cient of 10 -3 ) was considered where the reduction in artificial dissipation with fluctu- 
ation splitting results in a significant accuracy improvement over DMFDSFV. Con- 
vergence times were compared, showing a speedup of 2.5 for fluctuation splitting over 
DMFDSFV on identical grids, using a point Gauss-Seidel relaxation. However, a grid 
convergence study shows fluctuation splitting has better resolution of the solution on 
a coarser mesh than DMFDSFV does on finer meshes, resulting is a speedup of 29 
for fluctuation splitting over DMFDSFV. 



Chapter 4 

Scalar Mesh Adaption 


Mesh adaption strategies, for use in conjunction with the discretization and solution 
procedures discussed in chapter 3, are sought for two-dimensional scalar problems. 
Mesh adaption seeks to reduce to a minimum the number of nodes needed to obtain a 
desired level of solution accuracy. Implicit in this objective is the optimal placement 
of the retained nodes for a “most-accurate” solution. 

Many authors[61, 62, 65, 66, 67, 94] prefer to maintain isotropic (or nearly so) 
elements while performing unstructured mesh adaption. Their belief is that the ac- 
curacy of the solvers degrades on anisotropic elements. However, the use of isotropic 
elements will lead to excessive mesh densities when resolving phenomenon with dis- 
parate length scales in multiple dimensions, e.g. plane shocks or boundary layers. 
A contrasting approach is to employ anisotropically stretched elements, a technique 
long used for structured meshes[57, 70]. More recently, highly anisotropic elements 
have been demonstrated with unstructured grids[71, 74, 91, 95]. 

Two primary techniques are commonly used to perform mesh adaptions. One is 
to perform a global remeshing [96], where the adapted mesh is unrelated to the initial 
grid. This method works well when the adapted mesh will be significantly different in 
character from the initial mesh and usually employs an efficient grid-generation code 
to accomplish the adaption, with the initial solution driving the placement of source 
terms for the grid generator. 
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An alternative strategy is to form the adapted mesh from a series of local opera- 
tions on the initial mesh, based upon the local solution. This approach is very efficient 
when only a smaller fraction of the total grid requires refinement. Also, local opera- 
tions can allow for more control over the adapted grid and can be more robust than 
global grid generators. On the other hand, global remeshings can produce smoother 
grid distributions than can local operations. The present dissertation will pursue the 
local adaption strategy, as opposed to global remeshings. As before, unstructured 
triangulated domains are considered with the solution values stored at the nodes. 

This chapter begins by introducing the current state-of-the-art for anisotropic local 
refinement for finite volume solution methods as advocated by Habashi[71, 74, 75, 76]. 
Then a consideration of the behavior of fluctuation splitting methods on general 
meshes is presented. From this, a framework for local, anisotropic mesh refine- 
ment is developed in conjunction with the multi-dimensional fluctuation splitting 
scheme. The behavior of the new fluctuation splitting adaption strategy has a feature- 
alignment flavor, rather than feature-clustering like current methods, implying a re- 
duction in the required number of nodes, and hence solution times, for non-linear 
problems. 


4.1 Curvature Clustering 

Current state-of-the-art for unstructured mesh adaption is based on a posteriori error 
estimation [61, 62, 65, 66, 67, 69, 71, 74, 75, 76, 97, 98, 99]. The error estimates 
can be derived either by looking at the leading truncation error terms in the spatial 
discretization or by considering the solution interpolation error. In practice, either 
approach reduces, for second-order-accurate spatial discretizations, to a check on the 
curvature of the solution. Nodal positioning is then driven to equate the scaled 
second derivatives along all edges. If isotropic cells are desired, the magnitude of 
the local curvature inversely dictates the element sizes, while for anisotropic adaption 
directional derivatives can be used to stretch elements. 
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Edge in bold 
is not candi- 
date for swap- 
ping. 


Edge in bold Edge has been 

is candidate swapped, 

for swapping. 


Figure 4.1: Edge swapping schematic. 


Habashi[76] defines the a posteriori error estimate on an edge to vary like, 

\E r \ 

or as “the edge length squared times the second derivative of the solution.” In the 
finite volume context, one simple and efficient way to construct an edge-based error 
estimate along edge 01 is, 



The cdgc-bascd finite volume scheme detailed in chapter 3 already performs the gra- 
dient computations, making the error estimate a trivial computation. 

4.1.1 Adaption Mechanics 

Edge swapping 

One simple method to improve a mesh is to swap edges between nodes, altering the 
local connectivity. If the triangles to either side of an edge form a convex quadrilateral 
(Figure 4.1) then that edge is a candidate to be swapped. If the error estimate for 
the swapped connectivity is smaller than for the current edge connectivity, then the 


d l U 


dr 2 


(4.1) 
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Figure 4.2: Initial, poorly aligned mesh and solution. Solution contours vary on (0,1) 
with 0.1 increments. 


edge is swapped. In practice, an error threshold is employed to avoid limit cycles in 
smooth regions of the solution. 

The edge swapping procedure is demonstrated in Figures 4. 2-4.4. The problem 
statement is linear advection of a shear discontinuity with A = (1, 1) and a uniform 
inflow on the left boundary U(—l,y) = 0 and uniform inflow' on the bottom bound- 
ary U(x,0) = 1. The DMFDSFV scheme is used and reconstruction is performed 
using the compressive Van Albada limiter. The initial mesh is poorly aligned, 
Figure 4.2(left), and the resulting solution shows significant contour spreading, Fig- 
ure 4.2(right). Successive global edge swapping sweeps, Figures 4.3 and 4.4, improve 
the solution, reducing the contour spread. 

Excessive dissipation is created in the bottom left inflow' cell for this case, where 
the diagonal is not swapped with the present routine. Node insertion, to be discussed 
presently, will be able to handle cells such as these, wfliere the edge swapping logic 
breaks down. 








76 


CHAPTER 4. SCALAR MESH ADAPTION 




Figure 4.6: Four edges are the minimum that must be connected to node 0. 


Point deletion 

If all edges connected to a node have an error estimate, Eqn. 4.2, less than a threshold 
value, then that node can be removed. Once a node is flagged for deletion, the number 
of edges emanating from that node are reduced via edge swapping, as depicted in 
Figure 4.5. Usually, the number of edges connected at a node can be reduced to 
three, though Figure 4.6 shows a case where the minimum number of connected 
edges at the node is four. 

The removal of one node also eliminates two triangular elements and three edges. 
A schematic of the actual removal process of an interior node is depicted in Figure 4.7 
for both the three and four edge connectivity cases. The removal of a boundary node 
with three edges connected, depicted in Figure 4.8, eliminates one triangular element 
and only two edges. 

The objective of point deletion is to minimize the number of grid points while 
retaining the same accuracy. In practice, bookkeeping of node/cell/edge numbering 
is the most difficult part of the point deletion process. 
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-A 



(a) Three edges connected. 


(b) Four edges connected. 


Figure 4.7: Schematic of node deletion. 




Figure 4.8: Schematic of boundary node deletion. 


A numerical example of point deletion is presented for the 45° linear advection 
example. Figure 4.9 shows a grid and converged DMFDSFV solution to which the 
edge swapping routine has been applied. The point deletion routine is then applied, 
resulting in the grid of Figure 4.10. The solution on the reduced mesh remains 
unchanged from the finer-mesh solution. 

Point insertion 

Point insertion can be considered to be the opposite procedure of the four-edge con- 
nectivity point deletion, shown in Figure 4.7(b). If the error estimate along an edge 
exceeds a given threshold, then that edge is flagged for refinement. A new node is 
added at the midpoint of the edge, two new cells are defined, and three new r edges 
are created. The aim of point insertion is to increase solution accuracy, albeit at an 
increased computational cost due to the additional mesh nodes. Use of the local error 
estimator allows for selective refinement and results in fewer nodes than would be 
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Figure 4.11: After point insertion. Solution contours vary on (0,1) with 0.1 incre- 
ments. 

produced by a global refinement. 

The familiar numerical example of this section is continued by applying the point 
insertion algorithm to the solution and mesh last shown in Figure 4.10. The locally 
refined mesh, with most of the refinement occurring near the inflow boundary where 
the greatest artificial diffusion is produced, is shown in Figure 4.11, where 13 new 
nodes have been introduced. The solution shows improved accuracy, particularly in 
the reduced spreading in the lower left corner. 

Nodal displacement 

The fourth local adaption technique considered is nodal displacement. Neither the 
number of nodes nor their connectivity are altered, but their spatial location is op- 
timized. The nodal displacement concept mimics the Gauss collocation quadrature 
points familiar in calculus, and hence seeks to improve solution accuracy without 
increasing computational cost. 

The approach taken here uses the spring analogy, similar to the techniques pre- 
sented by Gnoffo[57] and Ait-Ali-Yahia[70] et al. The springs are taken to be the 
mesh edges. The spring constants are the edge error estimates. A minimization is 
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Figure 4.12: After nodal displacements. Solution contours vary on (0,1) with 0.1 
increments. 


then sought for a potential energy formed as, 

/*" 2 I ^ T ^ ^oj j is all distance 1 neighbors to node 0 (4.3) 

j 

The minimum is obtained by setting the derivative to zero, giving the requirement, 


J2\ E r\j f 0'j =0 (4.4) 

j 

where O' indicates the new, minimum energy position of node 0. Holding the edge 
error estimate constant during the displacement, an updated position vector can be 
obtained, 


{p)'j + ftyo) — 0 

j 


r jo' 


^hj\^r\j r 0j 


(4.5) 

(4.6) 


Equation 4.6 gives the position vector pointing from the current location of node 0 
to its adapted location. This adapted location is an attempt to optimally equate the 
scaled error estimates over all edges connected to the current node. Since the error 
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estimates are held constant during the local adaption step, an iterative process is 
required to obtain the optimal nodal positions. 

The nodal displacement routine is applied to the results in Figure 4.11, completing 
one global iteration of the full mesh adaption suite for unstructured meshes, consist- 
ing of nodal displacements, insertions, and deletions and edge swapping. Solution 
improvement, shown in the tighter contour clustering in Figure 4.12, particularly 
toward the upper right, is demonstrated using the point movement technique. 

4.1.2 Adaption Procedure 

Following the recommendations of Dompierre[71] et al. and Habashi[76] et al, the 
nodal displacement technique is used as a smoother between applications of the other 
three adaption operations. Also, close coupling between the solver and the mesh 
adaptor is to be maintained. 

The overall solution strategy therefore follows: 


1. Solve on initial mesh. 

2. Swap diagonals and iterate solver. 

3. Move nodes and iterate solver. 

4. Insert nodes and iterate solver. 

5. Move nodes and iterate solver. 

6. Delete nodes and iterate solver. 

7. Swap diagonals and iterate solver. 

8. Move nodes and iterate solver. 


Steps 4-8 are repeated until convergence of the entire process. 

The entire solution-adaptive procedure is applied to the 45° shear problem using 
the DAfFDSFV solver with the Van Albada limiter. The initial grid, containing 
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Figure 4.14: Adapted mesh with greatly improved solution. Solution contours vary 
on (0,1) with 0.1 increments. 
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121 nodes, and the solution on this mesh are shown in Figure 4.13. After only four 
complete adaption cycles, the mesh size has been reduced to just 103 nodes while 
capturing a very accurate solution, as seen in Figure 4.14. While the mesh adaption 
is highly anisotropic, the solution remains well behaved. 

A demonstration of the method is also performed for a case with a variable ad- 
vection velocity, namely the circular advection problem. The inflow profile contains 
a linear ramp between U = 0 on .t < —0.7 and U = 1 on x > — 0.6. The initial mesh 
is the same as for the shear case, containing 121 nodes. The solution on this mesh 
is excessively smeared, Figure 4.15. Three mesh adaption cycles were performed, re- 
sulting in a final mesh having 146 nodes. A significant improvement is again seen in 
the resulting solution, Figure 4.16. 


4.2 Characteristic Alignment 

4.2.1 Exact Advection Solution 

While anisotropic adaption to reduce a posteriori error estimates represents current 
state-of-the-art for unstructured-grid finite volume algorithms, it was shown in chap- 
ter 3 that there exists a better scheme for multi-dimensional advection than DMFDS- 
FV, namely fluctuation splitting. The question arises as to whether an edge-based 
error estimate, which is by the definition of an edge a locally one-dimensional indi- 
cator, is an appropriate adaption criteria for a truly multi-dimensional scheme. This 
section explores the development of candidate multi-dimensional adaption criteria 
that fit within the framework of the four local anisotropic adaption operations — edge 
swapping, point insertion, point deletion, and nodal displacement. 

As a starting point, posit that an optimal mesh is one with the least number of 
nodes on which a positive scheme can obtain a solution free from artificial dissipation. 
The artificial dissipation terms for fluctuation splitting, (f> x and <j>'" in Eqn. 3.29, 
were seen to scale with the limited fluctuations, (f>* X and cj>*" in Eqns. 3.26 and 3.27, 
respectively. For a locally positive scheme the limited fluctuations are bounded by 
the total cell fluctuation. Thus, driving the cell fluctuation to zero will eliminate 
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Figure 4.16: Converged DMFDSFV solution and mesh after three mesh adaption 
cycles, A = (y,—x). Solution contours vary on (0,1) with 0.1 increments. 
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artificial dissipation production, and vice versa. 

Combining Eqns. 3.22, 3.24, and 3.25, the cell fluctuation is, 

cb = -|nx -A(U 2 -U x ) + ~h ?J -A(U ?> - U 2 ) (4.7) 

or, 

4> = \ [hniUi - (£ih x + £ 3 h 3 )U 2 + i 3 h 3 U 3 ] A (4.8) 

Noting the vector identity ^) 3 =1 = 0’ the fluctuation can be expressed in compact 

notation, 

1 3 

9 =-Y J U i i i hi-A (4.9) 

;=i 

Recall from the notation for an element’s geometry, Figure 3.3, that edge i is the edge 
on the opposite side of the triangle from node i. If the cell-averaged flux Jacobian, ,4, 
is parallel to side i, then hi- A = 0. In the context of mesh adaption, assume that 
an edge, say without loss of generality edge 3, is aligned to be parallel to A, then 
from Eqn. 4.7 it is clear that the fluctuation, and hence the artificial dissipation, goes 
to zero as U 2 —> U x . Note that this condition for a vanishing fluctuation holds true 
irrespective of the value U 3 . 

The known solution to steady advection is that the solution remains invariant 
along the characteristics. A coupled solution-adaptive strategy is apparent — the adap- 
tion aligns one edge of each cell with the discrete counterpart to the characteristic, 
namely A, and the fluctuation splitting scheme converges to an exact advection so- 
lution by equating each downstream node of an aligned edge with the value of the 
upstream node. Efforts along this tack follow, though not without some lingering 
questions. What of the non-linear problem, where characteristics vanish at a shock? 
How can systems of equations be handled with multiple or imaginary characteristics 
defined on each cell? What about combined advection-diffusion problems where the 
exact solution no longer remains invariant along characteristics? 

By way of demonstration that in fact fluctuation splitting does compute the exact 
solution to multi-dimensional linear advection problems when the grid is aligned such 
that each triangle with non-zero gradient has one edge parallel to the cell-average 
advection velocity, the two cases employed with the curvature clustering adaption are 
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considered. The 45° shear case can be captured using a mesh with only 6 nodes, none 
of which are on the interior. The optimal mesh and fluctuation splitting solution 
are shown in Figure 4.17. Contrast these results with the mesh-adapted DMFDSFV 
solution to this problem in Figure 4.14. Fluctuation splitting in this case produces 
a more accurate, diffusion-free solution on a mesh 17 times smaller. A DMFDSFV 
solution is performed on this optimal fluctuation splitting grid, yielding the results 
of Figure 4.18, where 30 percent of the shear discontinuity has been smeared by the 
outflow boundary. 

For the circular advection problem a mesh is constructed containing 10 nodes, an 
order of magnitude fewer than were in the adapted DMFDSFV results of Figure 4.16. 
The fluctuation splitting mesh and solution are shown in Figure 4.19, where the inflow 
profile is exactly mapped to the outflow boundary. The mesh is called ‘optimal’ 
for this case because it contains the fewest points required to convect a diffusion- 
free solution. However, the solution representation on the interior of the domain is 
exact only in a discrete sense. Note that while the DMFDSFV solutions used the 
compressive Van Albada limiter, the fluctuation splitting solutions were performed 
with the non-compressive Minmocl limiter. The DMFDSFV solution on the optimal 
fluctuation splitting grid is shown in Figure 4.20, where it is clear that significantly 
more nodes will be required to achieve reasonable accuracy. 

These test cases suggest a tremendous potential for highly accurate fluctuation 
splitting solutions on extremely coarse meshes, when properly aligned. 


4.2.2 Adaption Strategies for Fluctuation Splitting 

The fundamental strategy for optimizing a mesh for fluctuation splitting solutions to 
linear-advection problems is to align one edge of each triangle with the cell-average 
characteristic direction. Translating this requirement into a systematic logic for local 
adaption is not straight-forward, particularly for a variable advection velocity when 
one edge cannot simultaneously be aligned with the characteristics in the elements to 
both the right and left of the edge. 
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Figure 4.17: Optimal mesh and diffusion-free solution using fluctuation splitting, 
A = (1, 1). Solution contours vary on (0,1) with 0.1 increments. 



Figure 4.18: DMFDSFV solution on optimal fluctuation splitting mesh, A = (1,1). 
Solution contours vary on (0,1) with 0.1 increments. 
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Figure 4.19: Fluctuation splitting optimal grid and diffusion-free solution, 

A = (y, — x ). Solution contours vary on (0,1) with 0.1 increments. 



Figure 4.20: DMFDSFV solution on optimal fluctuation splitting mesh, A = (y, —x). 
Solution contours vary on (0,1) with 0.1 increments. 
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(a) Misaligned mesh. 


(b) Aligned mesh. 


Figure 4.21: Edge swap options for linear advection with fluctuation splitting. 


Edge swapping 

Consider the 45° shear advection on a unit square, Figure 4.21. For this problem 
A = A = (1, 1) for all cells. The choice of diagonal perpendicular to A, Figure 4.21(a), 
constitutes a poorly aligned grid. The other choice of diagonal, Figure 4.21(b), creates 
a perfectly aligned grid. For the misaligned grid the cell fluctuations can be computed 
using Eqn. 4.9, 

d>A = \{2U X - U-2 - U 4 ) <Pb = \(U 2 + U 4 - 2 U 3 ) (4.10) 

and their sum is 4>a + (f>B = U 4 — C/ 3 . 

For the grid-aligned case, the fluctuations are, 

<ba = |(C/r - U 3 ) = <p b = ( ^hi (4.11) 

which leads to the observation, 

4*A T (pB = (ba + 4>b = 2^> a (4.12) 

Clearly, swapping the diagonal does not in itself reduce the total fluctuation. Yet it 
is known that fluctuation splitting will obtain the exact solution on the aligned grid 
and a diffused solution on the misaligned grid (recall the results of section 3.2.4). 
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Still seeking a link between mesh alignment and fluctuations, consider the RMS 
of the fluctuations. It is postulated that the RMS of the fluctuations is less on the 
aligned mesh than on the misaligned mesh. The proof proceeds from Eqn. 4.12, 


04 + 0B 
(j) 2 A + 20 a4>b + 0b 
04 + 0B 


= 2 4>a 

= 40* 

= 4 (j) 2 a - 2 <j) A (f>B 

— 4 0^0b) + 20, 
04 + 0b x 


= \ 


1 

2 


— 44>a<Pb 


+ 20 * 


(4.13) 

(4.14) 

(4.15) 

(4.16) 

(4.17) 


04 + 0B — ^(04 — 20A0B + 0b) + 20* (4.18) 

= ^(0.4 - 0b) 2 + 20^ (4.19) 

>202 = 02 + 02 (4.20) 

The use of the RMS of fluctuations to guide adaption for scalar problems has 
been proposed by Roe[32, 77]. The fluctuation splitting solver Roe employs is of the 
non-positive type, and for general solutions suffers from extreme dispersion errors. 
Also, nodal displacements are the only adaption operation demonstrated. 

A demonstration of the present edge swapping is performed for the 45° shear 
problem. The initial grid and fluctuation splitting solution are shown in Figure 4.22. 
The fluctuation splitting solution on this mesh is similar to the one obtained using 
DMFDSFV in Figure 4.13. However, after just a single edge swapping sweep the 
exact solution is obtained by fluctuation splitting, Figure 4.23. DMFDSFV was not 
able to match this accuracy from fluctuation splitting even on the highly adapted 
mesh from Figure 4.14. 


Point deletion 

Point deletion for a fluctuation splitting grid mimics the procedure previously dis- 
cussed in the finite volume context. If all elements connected at a given node have 
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Figure 4.22: Starting mesh and converged fluctuation splitting solution, A = (1.1). 
Solution contours vary on (0,1) with 0.1 increments. 



Figure 4.23: Exact solution obtained by fluctuation splitting after one edge swap 
cycle, A = (1,1). Solution contours vary on (0,1) with 0.1 increments. 
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small or vanishing fluctuations, then that node is a candidate for deletion. Edge 
swapping is performed to reduce the number of edges, and hence cells, connected to 
the node, under the restraint that the RMS fluctuations of the swapped configura- 
tions must remain small. If the number of cells connected at the node can be reduced 
to three, then the node can be deleted if the fluctuation remains small in the triangle 
formed by the agglomeration of the three connected elements. If four elements remain 
connected at the node, the point may be deleted if the RMS fluctuations remain small 
in either of the two possible agglomerations from four triangles to two triangles. 

When the point deletion algorithm is applied to the advected shear results seen 
in the edge swapping demonstration, Figure 4.23, the minimum optimal grid and 
associated exact solution, presented in Figure 4.17, are obtained after a single loop 
over the domain. 

Nodal displacement 

A scheme for displacing nodal coordinates to minimize RMS of fluctuations about 
the node has been presented by Roe[32, 77]. The development presented by Roe 
uses the same minimization scheme to evolve the solution as is used to drive the 
nodal positioning. That type of fluctuation splitting solver has a central difference 
flavor and is not monotonic. The present procedure incorporates the upwind, non- 
linear fluctuation splitting algorithm detailed in section 3.2.1 with some of the mesh 
movement strategies of Roe. The development also bears some resemblance to the 
adjoint equation adaption schemes of Anderson[98] and Venditti[99]. However, the 
adjoint schemes seek to evaluate changes in global parameters with respect to an 
error estimate computed on a coarsened mesh for a finite volume solver, while the 
present technique is entirely local and does not employ a coarsened sub-mesh and is 
specifically for a fluctuation splitting discretization. 

At a given node, the nodal displacements are computed as a first step and then 
the solution is updated at the new nodal location via a local point-implicit inversion. 
In this manner a global mesh movement sweep can be accomplished in conjunction 
with a single Gauss-Seidel iteration of the solver. In this section the current node to 
be moved is globally numbered node i. Within each triangular element the nodes are 
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locally numbered 1-3. Derivatives in y are omitted when they exactly follow from the 
derivatives in x. 

The functional to be minimized is defined at the node as a sum of contributions 
from all cells surrounding the node (equivalent to Eqn. 13, p. 248 in [77]), 

Ti = 2 ^ T “ T ^ T = 2 X] (4-21) 

VT "" T 

The weighting factor, Sj, is a positive scalar, generalizing to a symmetric, positive 
definite matrix for systems. The functional is thus a sum of positive semi-definite 
contributions from triangles containing the current node. 

The derivative of T with respect to a nodal coordinate is, 


dT 

dx 


=E 


(pj d^T , m , d(f>T 


2 dx, 


+ 


dx, 


(4.22) 


Note that the derivatives in Eqn. 4.22 represent the change in solution values as 
the discrete mesh is perturbed, and as such are to be interpreted in the context 
of variational calculus, and not as spatial gradients according to the more-familiar 
multi-variable calculus. 

The minimization of T can be performed using a fixed-point iteration to force the 
derivatives to zero, 

or 

dxi 

dr 

dx, 


X n+l = r" 


A n x; =- — 


(4.23) 

(4.24) 


Equation 4.24 can be combined with Eqn. 4.22 in the form of a distribution method 
of steepest descent, 


A”,. (4.25) 

l dxi dxi 

Convergence can be enhanced over the fixed-point iteration by using a Newton 
scheme. Expanding the gradient in an approximate Taylor series, 
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leading to the form, 
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In [32], Roe suggests neglecting the off-diagonal terms in Eqn. 4.28. 
Second derivatives of the objective function are, 
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(4.30) 


The fluctuation can be manipulated from Eqn. 4.9 as, 

<t> =\ [Ui(y 3 - 2 / 2 , x 2 - £ 3 ) + £' 2 ( 2/1 - 2 / 3 , x 3 - ari) + U 3 {y 2 - yi, - a: 2 )]-l 
= l [(U 2 ~ U 3 )(y 1 , -aq) + (££3 - ££0(2/2, -* 2 ) + (Efi - U 2 )(y 3 , -x 3 )]-X (4.31) 


The spatial derivatives take the form, 
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(4.32) 

(4.33) 


The derivatives of the cell-averaged flux Jacobian will depend upon the particular 
flux function. However, since A is a weighted average of three nodal values, 
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p. on (4.34) 
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For uniform advection, ^ = 0. For circular advection, J4- = ( 0 , — |) and = (|, 0). 
For non-linear problems, in general, 
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In keeping with the Gauss-Seidel update philosophy, only the derivative of the 
solution value at the current node is retained in the last term of Eqns. 4.32 and 4.33. 
That is, if the current node is designated node 2 of the triangle under consideration, 
then is retained while ~ W 2 - ~ 0 is assumed. 

Evaluating directly from the high-resolution non-linear fluctuation splitting 
scheme is impractical. Limiters such as Minmod are not continuously differentiable, 
and while the Van Albada limiter is differentiable its use does not lead to a con- 
venient explicit form from which to evaluate As an approximation to for 
the non-linear scheme, derivatives are sought using linear distribution schemes. Two 
linear choices for fluctuation splitting are to use a linearity preserving (second-order 
spatial accuracy), non-monotonic distribution or an upwind, monotonic first-order 
distribution. The linearity-preserving, non-monotonic scheme is of the Lax-Wendroff 
type[100], and tends to produce dispersion waves in response to nodal displacements. 
This behavior tends to under-predict the change in solution value at the perturbed 
node relative to the fluctuation splitting scheme employed herein. 

The linear upwind scheme is obtained from the present fluctuation splitting scheme 
by discarding the limiter. This scheme exhibits a dependency on the numbering of 
nodes within a triangle. To alleviate this dependency, the approximation to is 
built by looping over all cells connected to node i and renumbering the nodes within 
each triangle so that the current node is designated as node 2 of the triangle. It 
is emphasized that the linear upwind scheme is not used in the calculation of the 
solution, but only employed to provide an estimate of the solution variation with 
respect to nodal displacements, providing the forcing functions for the mesh adaption. 

The distribution to node 2 of a triangle using the linear upwind scheme is obtained 
from Eqns. 3.22-3.30, with ft^ = ft and ft" = ft, as, 

a + (U \ - U 2 ) - P~(U 3 - U 2 ) = a + L\ + (/T - a + )U 2 - (HI ft (4.36) 


where, 


a = 


a±|al 1 ± sign(a) 

— - — = a — ■ — 


0* = (4.37) 


Assembling all contributions from the surrounding cells of the form in Eqn. 4.36 and 
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solving for the steady state value of the current node yields, with U 2 = U t . 




(4.38) 


The variation of the nodal solution with respect to nodal displacement can now be 
evaluated as, 


®Uj + _ o-\ _ v - [tj ^ a+ _ jt 

Rrr- ^ _ ( Ul Fir. 3 r)v ■ 


Ui'E 


da+ dfj~ 
dx; dxi 


(4.39) 


Recall that the solution at surrounding nodes is held fixed during displacements of 
the current node. 

The derivatives of the functions in Eqn. 4.37 are defined, 


0 a < 0 
^ a>0 

OX — 


dij- f 0 < 0 


0 (5> 0 


(4.40) 


Further, Eqn. 3.22 leads to, 


1 „ „ dA 

— + -l|/q •- — 
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i y 1 dA 
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A* 1 , dA 

— — + -IlTli-— — 
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,4 X 1 , „ dA 

~--C nr di 


(4.41) 


(4.42) 


For non-linear problems an implicit relation arises for and One option is to 
neglect the derivatives of A in Eqns. 4.41 and 4.42. Another option would be to lag 
these same terms from the previous iteration level. 

Second derivatives of the fluctuation with respect to variation of a nodal location 
follow from Eqns. 4.32 and 4.33, again incorporating the approximation ~ ~0, 


, . rr.dAy _ dA dU 2 i 2 A ?d 2 U 2 vZ TJ L „ d 2 A 

(tl “ m + <2n2 'a^ + l"' 2 ' A ~dx( + T, L i2 n n S4 


r TT r-,dA X „ 8AdU 2 h . ~x8 2 U 2 U „ ^ A 

(U 3 — lJ \ ) — h i 2 n 2 ■ — — ~T b —772 -^02 + bj —rij ■ — J 

dy 2 dy 2 dy 2 2 8y{ J 2 8y{ 


(4.43) 


(4.44) 
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d 2 (f> 
dx 2 dy 2 




dA dlh dA dU> 
+ dx 2 ()y 2 



d' 2 A 
dx 2 dy 2 


(4.45) 


Second derivatives of A are developed in an analogous manner to its first derivatives. 

The second derivatives of the solution at the current node with respect to varia- 
tions in the position of the current node can be developed from the first derivative 
expression in Eqn. 4.39, again with the approximation ~ ~ 0, 


d 2 U t 
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(4.46) 


d£_\ 

dyi ) 

) (4-47) 


The remaining derivatives to be specified follow from Eqns. 4.41 and 4.42, 


d 2 a 
dx 2 


d 2 /5 
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For E, Roe[77] chooses E T = A-, for which the derivatives are, 

9B t _ 

dx,, Sj dx, 

d 2 Ej 2 fdS T \ 2 1 d 2 Sj 

dx 2 Sj \ dx, J Sj dx 2 
d 2 S T _ 2_dSrdST 1 d 2 S T 

dxidyi S\ dxi dy, S'j dx t dy , 


(4.50) 

(4.51) 

(4.52) 


The area of a triangle is, 

Sj = \ M?/2 - ys) + x 2 (y 3 ~ yi) + xz{yi - y 2 )] 

= bi(^2 - x 3 ) + y 2 (x 3 - xi) + y 3 {x 1 - x 2 )\ 


(4.53) 


leading to the derivatives, 

dSj 
dx 2 

and, 


17.3 ~ .Vi 
2 


dSj x\ — x 3 
dy 2 2 


(4.54) 


xs 1 = xs 1 = yy^ = 0 

dx 2 dy 2 dx 2 dy 2 

This choice of weighting emphasizes fluctuations on the smaller cells. 

An alternative is to weight all cells equally, with S T =1. Two other obvious choices 
for weighting are Et=5t, emphasizing fluctuations on the larger cells, and ER = 
for even stronger emphasis on the smaller cells. The derivatives for this latter case 
are, with Eqn. 4.55, 


d^L x 2 dSj 

dxi Sj dxi 

d 2 E T _ _ 6 _ fdSj \ 2 

dx 2 S'! \ dx, J 

d 2 E T _ 6 dS T dS T 
dx^ dxi dy. 


(4.56) 

(4.57) 

(4.58) 


The fluctuation splitting node movement techniques are evaluated for the familiar 
circular advection problem of the present chapter. One cycle of point deletion results 
in the grid and solution of Figure 4.24. The grid contains 70 nodes. 
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Figure 4.24: Starting mesh and solution for demonstrating fluctuation splitting node 
movement schemes. Contours vary on (0,1), with 0.1 increment. 



Figure 4.25: Mesh and solution after four fixed-point mesh adaptions with fluctuation 
splitting. Contours vary on (0,1), with 0.1 increment. 
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Figure 4.26: Mesh and solution after four diagonal Newton mesh adaptions with 
fluctuation splitting. Contours vary on (0,1), with 0.1 increment. 

Applying the fixed-point iteration strategy for four node movement cycles results 
in the mesh and solution of Figure 4.25. This fluctuation splitting solution shows 
a somewhat improved accuracy over the optimally adapted DMFDSFV solution of 
Figure 4.16, despite a mesh with fewer than half the number of nodes. The fixed-point 
procedure is found to be a slowly converging technique. To improve the response, 
local iterations and/or over-relaxation can be applied. Unfortunately, the amount of 
over /under-relaxation for optimal convergence is problem and grid dependent. For 
the results of Figure 4.25, an over-relaxation factor of 10 was used with five local 
iterations per global adaption, with Hj=l. 

Employing the full Newton strategy to improve the convergence of the minimiza- 
tion procedure has not been encouraging for this case. Numerical stiffness due to a 
vanishing determinant in Eqn. 4.28 is the culprit, a condition that is likely to worsen 
for finer grids. 

However, the diagonalized Newton scheme proves robust and fast for this case, 
without the need for either local iterations or a tunable over/under-relaxation param- 
eter. Four cycles of the diagonal Newton scheme (S T =1) yield the grid and solution 
in Figure 4.26. The mesh adaption is clearly more aggressive than for the fixed- 
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point scheme, achieving slightly improved resolution of the solution at the outflow 
boundary. 


Point insertion 

Additional nodes are created by subdividing edges, in a manner similar to that de- 
scribed in section 4.1.1. However, instead of using an edge-based error estimate, the 
fluctuations in the cells to either side of the edge are used to trigger mesh refinement. 
If the magnitudes of the fluctuations in both cells, or the only cell for a boundary 
edge, exceed a threshold, then the edge is split by adding a new node at the edge 
midpoint. 

The various operations for performing mesh adaption in the fluctuation splitting 
context are combined as a series of sequential steps to form a complete adaption cy- 
cle. The steps are arranged in the same order as is enumerated in section 4.1.2. The 
complete fluctuation splitting mesh adaption is applied to the circular advection test 
problem for a single cycle. The initial mesh and solution were previously shown in 
Figure 4.15. The resulting mesh, containing 96 nodes, is shown in Figure 4.27(a). The 
associated solution, showing very good accuracy, is in Figure 4.27(b). The fluctuation 
splitting solution-adaptive procedure achieves slightly better solution resolution in a 
single adaption sweep than the state-of-the-art finite volume solution-adaptive pro- 
cedure was able to obtain in three adaption cycles on a mesh with 50 percent more 
nodes. Thus, the present fluctuation splitting adaptive scheme is more than three 
times as fast as the finite volume scheme which utilizes a posteriori error estimates. 

Note, however, that the ‘optimal’ fluctuation splitting mesh from Figure 4.19(a) is 
not achieved with only one adaption cycle. In fact, it appears that for this variable- 
velocity advection problem the local adaption strategy falls short of achieving the 
globally-optimal mesh. A retreat is therefore made from the goal of globally-optimal 
mesh adaption, and instead the present method seeks a consistently-improved solu- 
tion/mesh combination with each adaption cycle. 
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(a) (b) 

Figure 4.27: Mesh and solution for circular advection problem after a single mesh 
adaption cycle with fluctuation splitting. Contours vary on (0,1), with 0.1 increment. 

4.3 Non-linear Advection 

The two solution-adaptive schemes, a finite volume solver with adaption based on 
edge error estimates and a fluctuation splitting solver with adaption to minimize 
fluctuations, are applied to the two-dimensional non-linear advection problem of sec- 
tion 3.2.4, repeated here, 

U t + UU x + U y =0 

U{-l,y) = U{0,y) = 0 U{x, 0) = -2x - l 

To avoid the problem of a doubly-defined node at (—1, 0) and (0, 0), the inflow corners 
are set to U = 0, and thus the expansion at the inflow is linearly distributed between 
the corner points and the first node in from the corner, along the y-axis. 

The DMFDSFV scheme with edge error estimates is applied without modification. 
The initial mesh, Figure 4.28(a), contains 121 nodes. The DMFDSFV solution on this 
mesh, Figure 4.28(b), is very diffused and exhibits grid-induced asymmetries. Three 
full cycles of solution adaption are applied, nearly doubling the mesh density to 
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Figure 4.28: Initial mesh and DMFDSFV solution for non-linear advection case. 
Contours vary on (-1,1), with 0.1 increment. 


237 nodes, Figure 4.29(a). The solution on this mesh shows a dramatic improvement, 
Figure 4.29(b), for shock thickness, shock speed, point of coalescence, and preservation 
of extremum in smooth regions. Note, however, that there is some asymmetry between 
the expansion fans toward the top of the domain, that the shock is not entirely 
straight, and that the compression fan begins to coalesce into a shock at y = 0.4, 
instead of at y = —0.5, the correct location. 

The fluctuation splitting adaption scheme is also applied for three cycles to the 
same problem and initial grid. The solution on the initial mesh is similar in char- 
acter to the DMFDSFV solution in Figure 4.28(a). The adapted mesh, containing 
206 nodes, and corresponding solution are shown in Figure 4.30. The adapted fluctua- 
tion splitting solution, using 14 percent fewer nodes, exhibits slightly greater accuracy 
than the adapted DMFDSFV solution, particularly in the expansion fan symmetry 
and shock coalescence point. One feature that is better resolved in the DMFDSFV 
solution is the extremum between the compression and expansion fans on the low r er 
right-hand side. The fluctuation splitting shock is wider at the coalescence point but 
then has comparable crispness to the DMFDSFV result and is a little straighter wdth 
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Figure 4.29: Mesh and DMFDSFV solution after three adaption cycles, non-linear 
advection case. Contours vary on (-1,1), with 0.1 increment. 



Figure 4.30: Mesh and fluctuation splitting solution after three adaption cycles, non 
linear advection case. Contours vary on (-1,1), with 0.1 increment. 
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increasing y. The fluctuation splitting shock at the outflow is ever so slightly offset 
to the right from x = —0.5. the correct location. 

The DMFDSFV and fluctuation splitting meshes, Figures 4.29(a) and 4.30(a), 
have similar character in the expansion and compression fans. The fluctuation split- 
ting adaption does not cluster as many points to the shock, resulting in overall 14 per- 
cent fewer nodes. This fact that the fluctuation splitting adaption scheme can resolve 
crisp shocks without excessive clustering normal to the shock may have favorable im- 
plications when considering fluid dynamic problems with extremely strong shocks in 
the vicinity of more subtle, though still important, features, such as an entropy layer. 

While adapting to the non-linear problem, it became necessary to include an 
eigenvalue limiter in the fluctuation splitting algorithm to prevent expansion shocks. 
The one-dimensional eigenvalue limiting of Harden and Hymen (see section 2.2.2) is 
extended to multiple dimensions by searching for expansions in the £ and i] directions 
separately. Combining Eqns. 3.25-3.29, the artificial dissipation terms can be recast 
as, 


/ = -|a|A f t/ + sign(cv)A/*, — sign(, 6)M i , (4.59) 


|a| and \fi\ are then limited as in Eqn. 2.46 with the parameter e (Eqn. 2.47) defined 

as, 


e («) 

e (P) 
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- max 
2 

1 

- max 
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0 , 

0, 


£ihi‘(X- Ax), 
^ 3^'3 ’ (^2 ~ A ), 


Iifii ■ (A 2 — A) 
^3^3 ' {A A3) 


(4.60) 

(4.61) 


respectively. 

As a final note to this section, the fluctuation splitting solver is applied to the 
adapted mesh previously used with DMFDSFV, from Figure 4.29(a). The fluctuation 
splitting solution on this mesh is shown in Figure 4.31, as a demonstration that the 
fluctuation splitting solver still works on a mesh adapted to solution curvature. 
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Figure 4.31: Fluctuation splitting solution on curvature-clustering mesh. Contours 
vary on (-1,1), with 0.1 increment. 

4.4 Advection-Diffusion 

The two solution-adaptive schemes are applied to the Smith & Hutton advection-dif- 
f'usion problem of section 3.4, starting on the 1928-node Mesh C. The finite volume 
adaption with clustering based on edge error estimates is applied without modifica- 
tion, resulting in a 619-node mesh. The RMS difference between the solution profile 
at the outflow boundary for the adapted DMFDSFV result and the grid-converged 
solution profile from the previous chapter is listed in Table 4.1 along with the RMS 
difference for the DMFDSFV solution on the unadapted mesh C. Two cycles of mesh 
adaption have been performed, reducing the required number of nodes by two-thirds 
while improving the RMS outflow resolution by 28 percent. 

The incorporation of physical diffusion into the adaption strategy utilized with the 


DMFDSFV Fluctuation splitting 
adapted 0.0208 0.0063 

unadapted 0.0288 0.0068 


Table 4.1: RMS difference of solution outflow profile relative to grid-converged result. 
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fluctuation splitting scheme requires some modifications. A check on the magnitude 
of the diffusive fluctuations in the surrounding cells is added to the criteria for deleting 
a node. For edge swapping, the RMS of the diffusive distributions are added to the 
RMS advective fluctuations to determine the optimal connectivity. An edge is split, 
adding a node, if the diffusive distributions sent to the end nodes of the edge exceed 
a threshold value. 

For nodal displacements, the advective fluctuation on a cell, Eqn. 4.31, is modified 
to include the diffusive distribution sent to the current node, <$ defined in Eqn. 3.51, 
repeated here, 


*! = -7^'EVm-f‘i 

4St u 


(4.62) 
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where the assumption ~ ~ 0 has already been applied. The second derivatives 

of the diffusive distribution follow as, 
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The derivatives of the cell-averaged diffusion coefficient scale like, 

dp 1 dy 2 d 2 p 1 d 2 y 2 

dx 2 3 <9.x 2 dx 2 3 <Vr 2 


(4.68) 


The variation of the dependent variable at the current node with respect to the 
nodal position also needs to be modified to account for the diffusion terms. The 
steady-state value of the current node, Eqn. 4.38, now becomes, 
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and the variation, Eqn. 4.39, is replaced by, 
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The second variations of the dependent variable, Eqns. 4.46 and 4.47, become, 
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The RMS differences of the solution profiles at the outflow boundary relative to 
the grid-converged result are in Table 4.1. Although the fluctuation splitting solution 
on the unadapted 1928-node mesh is good, one cycle of mesh adaption reduces the 
number of nodes by a factor of 2.8 to 695 nodes, while still producing some (7 percent) 
improvement in accuracy. 


4.5 Recapitulation 

Current state-of-the-art for anisotropic unstructured mesh adaption based on a poste- 
riori error estimates is implemented in an edge-based structure in conjunction with a 
finite volume solver. This type of adaption results in meshes where the node densities 
are clustered to regions of high curvature in the solution. Significant improvement in 
solution accuracy is verified using this technique on scalar model problems. 

Recognizing the remarkable property of the discretized fluctuation splitting scheme 
that multi-dimensional advection can be solved exactly when one edge of each cell 



110 


CHAPTER 4. SCALAR MESH ADAPTION 


is aligned with the characteristic direction, a different mesh adaption scheme is pro- 
posed. Retaining the mechanics for performing only local operations, i.e. point inser- 
tion/deletion, edge swapping, and nodal displacement, a solution-predictive approach 
is chosen in favor of a posteriori curvature clustering. The concept of aligning cell 
edges with characteristic directions is generalized as a minimization procedure to al- 
low extension to diffusion problems and systems. Extending this process to non-linear 
problems leads to a multi-dimensional implementation of eigenvalue limiting. 

It is seen that while performing a series of local optimizations does lead to glob- 
ally improved solution accuracy and reduced grid sizes, in general a truly globally 
‘optimal’ mesh is not achieved in a reasonable number of adaption cycles. However, 
the solution-predictive adaption in conjunction with the fluctuation splitting scheme 
does provide moderately more accurate solutions on smaller meshes for comparable 
number of adaption cycles versus the DMFDSFV solver with adaption driven by error 
estimates. 

Considering extensions to three-dimensional hypersonic flow applications, perhaps 
the most promising difference between the two adaption strategies lies in their treat- 
ment of shocks. For the a posteriori adaption, the number of nodes clustered to the 
shock grows as the shock strength grows, which could lead to a bow shock dominating 
the adaption for hypersonic problems. In contrast, the minimization of fluctuations 
tends to merely align the grid with the shock, leaving the points outside the shock 
largely unaffected. 



Chapter 5 

Two-Dimensional Systems 


5.1 Overview 

The fluctuation splitting and DMFDSFV discretization schemes detailed in chapter 3 
for scalar advection-diffusion problems are extended in the present chapter to the 
system of equations governing the compressible fluid dynamics of a dilute, ideal gas 
for two-dimensional and axisymmetric flows. First, the mass, momentum, and energy 
equations for an inviscid fluid, termed the Euler[86] equations, are formulated in both 
the fluctuation splitting and finite volume contexts. Then, the system is extended to 
include the effects of viscosity and conductivity as the Navier[87]-Stokes[88] equations. 
The formulation of the upwind fluctuation splitting scheme for the two-dimensional 
and axisymmetric Navier-Stokes equations is a leading-edge research area. 

The chapter concludes with verification and validation of the schemes. A brief 
discussion of the coding strategy is included, immediately followed by a verification 
of the code using the test cases and methodology of Shinghal[101] and Roache[102]. 
The validation cases range from the incompressible flat-plate to a Mach-17 cylinder. 

5.2 Inviscid Formulations 

Spatial discretizations are developed for the Euler equations using both the DMFDS- 
FV and fluctuation splitting schemes. The two-dimensional development is extended 
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to cover axisymmetric problems through the additional appropriate radial derivatives, 
treated as source terms. 

The non-dimensional system of equations are taken from appendix B, Eqn. B.9, 
with F w = B'” = 0 as, 


m a U t + V-(tz7 0 F*) =wB'‘ 


(5.1) 
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(5.4) 

m is a logical switch between two-dimensional (m = 0) and axisymmetric (m = 1) 
equations and m a = 1 — m + my is toggled by m between 1 (m = 0) and y (m = 1). 


5.2.1 Finite Volume 


Integration of the dependent variables over the control volume about node i is per- 
formed as, 



m a \J t dil 


m a Sj\J 


(5.5) 


For two-dimensional W a = 1, while for axisymmetric m a can be either taken as 
= yi , for mass-lumping to the node, or as the y-value of the centroid of 
The axisymmetric source term, B', is simply evaluated at the node as, 



(5.6) 
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The discretization of the fluxes follows the development of section 3.2.1, with the 
only modifications being the inclusion of w a (equal to 1 for two-dimensional or the 
y - value of the quadrature point on the face for axisymmetric) and the definition of 
the artificial dissipation in Eqn. 3.8, which is, 

* =i|l-n|(U R -U L ) (5.7) 

where by convention the right state is to the outside of the control volume while the 
left state is to the inside. 

In an analogous manner to the conservative one-dimensional linearization of sec- 
tion 2.4.1, the parameter vector, Z = ^fp [1, u , v, H] T , is linearly averaged, 
Z = |(Z L + Z R ), to provide the quantities, 


Z-2 Z s ~ Z\ 

u = v = H = 

Z\ Z\ Z\ 

and the Roe-density remains from Eqn. 2.64, p = ^/PlPr. 

The projected flux Jacobian is decomposed as in Eqns. B.37-B.40, 


A-n| = X|A|X _1 


(5.9) 


with A given in Eqn. B.38 and X given in Eqn. B.40. 

The product X - 1 (Ur — U[_) is expressed analogously to Eqn. 2.68[23], 

2 a 2 dp — 2 dP 
2 a 2 (n x dv — n y du) 
dP + pa dV 
dP — pa dV 

where the projected velocity is V = V-n and the averaged speed of sound for a perfect 
gas is, 

5 2 =(7-l )(ff-^) (5.H) 
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5.2.2 Fluctuation Splitting 


In the fluctuation splitting context, the parameter vector is taken to vary linearly over 
each element. The assumed linear variation limits the scheme to second-order formal 
accuracy, but will allow for a positive formulation when combined with a non-linear 
distribution. Caraeni, Caraeni, and Fuchs[103], basing their work on the published 
formulations from the present report, use a quadratic variation of the parameter 
vector to derive a formally third-order scheme. That scheme, however, cannot be 
positive, a critical requirement for the robust computation of strong shocks. 

For a perfect gas, changes to the conserved variables can be related to changes in 
the parameter vector as, 


d\J =\J z dZ 

~‘2Z 1 0 0 0 " 

Z-2 Z\ 0 0 

u? - 

Z> y 0 Z\ 0 

_iz 4 \z x _ 

Integration of vo a \5 t over an element leads to a mass matrix, 


(5.12) 

(5.13) 


/ zu a U t dO = f WaU z Z t di} (5.14) 

Jn Jn 

If mass-lumping to the nodes is employed, introducing temporal, but not spatial, 
errors, Eqn. 5.14 can be distributed to the nodes as in Eqn. 3.14, so that the sum of 
all contributions to node i equals zu ai S t U H . 

While some authors insist on upwinding source terms[3, 104], the present analysis 
considers an upwind distribution to be inappropriate for the axisymmetric source 
terms, which arise from purely geometric manipulations. The axisymmetric source 
term can be distributed to the node in a mass-lumped analogy as, 


a, i ZU 


(5.15) 


which is equivalent to the finite volume treatment of Eqn. 5.6. A modification of this 
distribution is to send contributions weighted by the averaged values, 

vJaA U it ^w^Bj + COE 


(5.16) 
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Figure 5.1: Subdivision of triangular element into three quadrilateral integration 
areas. Dashed lines are the median-dual mesh. 


A more rigorous treatment integrates the source term analytically, based on a 
linear variation of the parameter vector. The only non-zero inviscid source term is, 


B' 





(5.17) 


The integration over the triangular element is divided into thirds along the median- 
dual boundaries, as in Figure 5.1, so that, 


Ft — IT T 1 1‘2 T 1 b 


(5.18) 


The subintegrals are then distributed to the nearest node. Notice that the subdivided 
integration elements, £ 2 i- 3 , are quadrilaterals, whereas the original element was a 
triangle. The distribution formula is thus, 

w ai Si\5 it i — zu I B 'dFli + COE (5.19) 

J 


The integration of the source term over is expanded in detail in appendix C.1.1. 
Integration of the inviscid flux is performed as, 
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The y-component of the flux function can be written in terms of the parameter vector 
as, 


z i 


Z1Z3 
Z2Z3 
- (ZiZ, 

Z-.\ Z,\ 



(5.21) 


A linear variation of the parameter vector over a. triangular element can be represented 
as, 

z(x,y) = 2s^ e ijk z j K x - x i)(vk - Vi) + (y- Vi)i x i - x u )] ( 5 . 22 ) 

where is the cyclic-permutation summation operator. The linear variation can 
also be written in the element-local (4, >]) coordinates, referring to Figure 3.3, as, 


z(£, n) = z L + f(z 2 - z ;[ k + Uz 3 - z 2 )y 


1 


1 


Z| + —\ZH + — -\Z // 

C3 tl 


(5.23) 


The domain is on 0 < i] < ^ and 0 < £ < £ 3 . The Cartesian coordinates map 
similarly, 


x^Al) = x i + ]~ \ x ^ + y\ x V ( 5 - 24 ) 

y(t,v) = yi + Y^y £ + v (5.25) 

t-3 tl 


Some general integration rules can be developed for linear variations over the trian- 
gular elements: 




— St 
= Sjx 


S T xy 


St 

4 




(5.26) 

(5.27) 


(5.28) 
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The cell-averaged value is, 

xi + x 2 + x 3 1 ^ 

x = g = 3 2^ x j (5-29) 

3 = 1 

The last term of Eqn. 5.20 is distributed to the nodes in a manner similar to the 
source term. 

<--w I F* y + COE (5.30) 

•to, 

The integration rules from Eqns. C.16-C.20 are used to evaluate the integral, leading 
to, 

f C _ _ _ 3 

/ ZiZ% dO t = — HZ 1 Z 3 + 11(Z[. Z 3 + ZiZ^) + 9Zi i Z^ t + Zj. (5.31) 
JQi L j = 1 

for the continuity equation. The integrals for the other governing equations follow 
directly from Eqn. 5.31. 

The remaining term to evaluate in Eqn. 5.20 is the inviscid fluctuation, 


- / w a V-FdO 

Jn 

— -FzA^Z — 4 ^ 3 'FzA^zj dQ 


where dF* = F^ dZ and, 


\ Z 2 Zi 0 0 j 



1_ 0 Z\ 0 Z 2 \ 

Z% 0 Z\ 0 

_ 0 ^3 ^2 0 
Z ^Zi -^Z 2 "^Z 3 ?=± 7 

/y ^ 7 7 7 J- 

0 0 ^4 Z3 


(5.33) 


(5.34) 
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The integration rule Eqn. 5.28 allows for the direct evaluation of Eqn. 5.32 as, 


1 — 1/— 1 -A vo a ^ \ 

^ = [f^ - i - 3 E ^F Z ,J j A { Z 


Noting that A = F z Zf, and, 


Z 77 — 


— H + (7 — 1) (« 1 2 + v 2 ) —2(7 — l)u —2(7 — l)v 27J 


with the tilde-averaged quantities defined as, 


U = U(Z) 


AfU = U/ AcZ 


A„U = U Z A„Z 


leads to the fluctuation expressed as, 


1 =, 1 z 1 > W a . - \ — — 

A-dA-.^A, A £ U 


1 ^ 1 / = 1 Cc7 a . - \ 

+ A-- A--^ A„U 


3 ' w s 

j=i 


= + (p V 


Equation 5.38 includes the approximation, 


A i - F A Z v 


A transformation of variables can be made to the auxiliary variables so that, 


4> — UV0 


(5.40) 
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with U w given in Eqn. B.61. The similarity transformations developed in section B.3 
lead to, 


0 2^a^lTl’ 


+ -^ahnr 


A 


\(X-\ v— A- 


i= i 


4 ' 


= — w. 


a (aA { W + i3A,w) 


j=i 


A f W 


A„W 


(5.41) 


0 s +0” 


with .4 = Wf/AUw and, 

a7w = w^u = w^a 5 z 


A„W = Wp A„U = A„Z (5.42) 


where Wj/ is given in Eqn. B.62 and. 


W z = Vp 


9 H _ 

7 h 7 h 


—U 
— V 


It 
7 h 


7 * 


7 


-^2 

-^3 

4^4 z 4 


1 

0 

0 


0 

1 

0 


7 U 

44 V 
7 

7-1 
7 - 


44 Z 2 

7 * 

7-1 

7a 2 ^ 

-44^1 

7 a z - 1 - 


0 

0 

0 

4 

0 

- 4 ^-Z 2 

7 

- 4 ^z 3 

4 ^Zj 

7 


using the relation, 


h = c p T = 

7-1 


The projection of A has the simple form, derived as Eqn. B.54, 


(5.43) 


(5.44) 



V 

0 

0 


n-A = 

0 

VI 

n T 

= 


0 

a 2 h 

V 



V 

0 

0 

0 


0 

V 

0 

a 2 n x 


0 0 

0 n x 
V n y 
a 2 n y V 


(5.45) 
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where the projected velocity is V = h-V. The generalized advection speeds are, 


« (A-lg^)] <*•«) 

0 = -§*„■ (5.47) 

incorporating the approximation, 

Aj =WfiAjiJ w (5.48) 

Linearity preservation for second-order spatial accuracy is obtained by limiting 
the fluctuations componentwise, 

= ¥j + 4>] <■(()/■ = (i - V’ ) (5.49) 

= 4>] - 4'{Qj) = # (1 - V> (Qj)) (5.50) 



with 


= diag ^1 - ip D ' 7 = diag(l - ^ (Qj)) (5.54) 

Upwinding is achieved through the introduction of artificial dissipation 

4 > = sign (a) 0 = — w a M a D f M“ 1 |o;| A ? W (5.55) 

$ " = sign(W " = -ro«Mf ! D’ , M^ 1 \j3\ A^W (5.56) 
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where M a = sign(a) and = sign(/3). The absolute values of the generalized 
advection speeds are developed using the following decomposition, which is exact for 
the two-dimensional equations but approximate for the axisymmetric form, 

a =^h 1 -A=jXAX~ 1 (5.57) 

where A and X are defined in Eqns. B.38 and B.57, respectively. The absolute value 
is then defined as, 

\a\ = l -±X\k\X~ l (5.58) 

Expressions for the sign of the generalized wavespeeds are developed from, 


|a| = M q q: 

M 0 = A'|A|A“ 1 Af“ 1 


(5.59) 

(5.60) 


leading to, 


Mq, = 


sign(V Q ) 

0 

0 

0 


M q = sign (V Q )I 

0 0 
n\ sign(V Q ) — n*nj[sign(V ( 

-n£n?sign(V a ) nfsign(V a ) 


if Vq. | > a 


an\ 


an\ 


nl 

a 

ni 

a 

o 


(5.61) 


if |V a | < a (5.62) 


where V a = h\-V . Similarly, defining |/3| as, 


101 = -^Ar|A|Af"‘ 


leads to. 


(5.63) 


M« = 


= sign (V#) 7 


sign(V^) 

0 

0 

0 


0 

n 3 sign(Va) 

-nWs sign(V^) 

av x . 


0 


if \%\ > a 


0 


-n^n|sign(V^) f 
nf F\gr\(V p) f 


ant 


(5.64) 


if |V^| < a (5.65) 
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where Vp = fi^-V. M Q and have the property, 

M - 1 = M 0 M^ 1 = M (i (5.66) 

Eigenvalue limiting for the suppression of expansion shocks can be introduced into 
the expression for |a|, Eqn. 5.58, by limiting |Ai| = |A 2 | = V|, |A 3 | = |V + a|, and 
| A 4 1 = \V — a\ in the manner of Eqns. 2.46 and 4.60. If the limited eigenvalue is 
expressed as, 

|A| Urn =|A| + Aa (5.67) 


then the additional artificial dissipation for eigenvalue limiting in the £ direction to 
be added to Eqn. 5.55 is, with A + = |(A A3 + A a J and A - = ^(Aa 3 — A A4 ), 



Aai 

0 

0 

0 

_ ix 

0 

rif A + + n\ A a . 2 

n\n{{ A+ - Aa 2 ) 

2 iA" 

a 

-w a - 

0 

n\n\{ A + - Aa 2 ) 

nf Aa 2 + n\ A + 

a 


0 

an* A - 

an\ A - 

A+ 

while the eigenvalue limiting in the t] direction takes the form, 


a Ai 

0 

0 

0 

h 

0 

nf A + + n\ Aa 2 

«3«s(A + - A.vJ 

Lk/\- 




2 &A - 

a 

2 

0 

n 3 n 3( A+ - A a . 2 ) 

t ?3 Aa 2 + n\ A + 


_ 0 

an .3 A - 

an\ A - 

A+ 


Ac W 


(5.68) 


A „W 


(5.69) 


The fluctuation is distributed to the nodes using a simple extension of the scalar 
distribution, Eqns. 3.30 and 3.31, 


SiU lt <- 
5*2 U 2i ■<— 

or in a more compact form, 


1«.* ! -/) + COE 
\(Y + A)+\ (A + <t>"') + COE 
\\P' -&") + COE 


(5.70) 


SiVi 


i( 3 - i )( + (-1 ytf) + (-4 + 5 ? - 


A COE 


i = 1,2,3 


(5.71) 
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where Eqn. 5.40 is used to define, 

0* 5 = il w <j>*\ </>'■ = \J W 0' (5.72) 

0*" =U w 4>*", 0'" = (5.73) 

At the 15 th AIAA Computational Fluid Dynamics Conference Prof. Philip Roe[105] 
suggested including the radial distance, w a , within the definition of the parameter 
vector. Further reflection reveals that sfw~ a is the proper term to consider. This 
approach has the potential to avoid some of the explicit integration by parts in this 
section and allow for cleaner expressions for many of the combined axisymmetric / two- 
dimensional relations, such as Eqn. 5.35. However, the flux function F is no longer 
just a function of state but also of position, so that the gradient can no longer be 
expanded as V-F = F^-VZ. which could lead to complicated expressions again. 
Also, the consistent incorporation of this redefinition of the parameter vector into 
the viscous discretizations is not straightforward due to the difficulty in expressing 
all terms as explicit quadratic products of the parameter vector. 


5.3 Viscous Formulations 

The complete set of governing equations for the motion of a two-dimensional or ax- 
isymmetric viscous fluid are presented in appendix B as Eqn. B.9, 

m a XJ t + V-(w 0 F l ) =V-{w a F v ) +wB ' 1 (5.74) 

Integrations of the conserved variables, inviscid flux, and inviscid axisymmetric source 
proceed as developed in section 5.2. 

Integrating the viscous flux over a triangular element leads to the viscous fluctu- 
ation (see Eqn. 3.47), 

<p v = I V-{w a F v )dn (5.75) 

Jn 

The nodal distributions are developed in a finite element sense by integrating by parts 
as the extension of Eqn. 3.49, 

0 " = <f v t zu 0 F v ■ fidT - I zu 0 F v ■ Vvi dil 

Jr J n 


(5.76) 
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For interior nodes the boundary integral in Eqn. 5.76 will sum to zero and the volume 
integral is integrated analytically for a linear variation of the parameter vector, 


VZ 


1 3 

' Zjijfij Yv{ 


j-i 


I) h % 
25V 


(5.77) 


0* = 


25 


T JQ 




The viscous fiux is, 


/ 




I! 


0 

r 


\kVT + Vt ) 


with the shear-stress tensor defined, 


r = /i 


v t f+ (v t v) t - |v-v/ 


(5.78) 


(5.79) 


(5.80) 


Struijs[34] et al. have shown that derivatives of primitive variables can be consis- 
tently represented in terms of the parameter vector gradients as, 


VV = Vz VZ 


(5.81) 


where, 


and, 


2 Z l 

0 

0 

0 

Zi 

C 

1 

Zl 

0 

0 

“# 

0 

1 

Z 1 

0 

zD-Z 4 

. 7 

- J ^-z 2 

7 

-^-Z z 

7 


V =V(Z) 


/ A 

12 

Zy 

la 

Vt [ziz.-lizi + zt)]) 


(5.82) 


(5.83) 
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Further, the consistent temperature gradient is, 


VT = 


VP PVp 


pR p 2 R Rp - \ 


( pVP — PVp 


(5.84) 


The viscous fluctuation is then distributed to the nodes as in Eqn. 3.52, 


<-0V + C7OE 


(5.85) 


An alternate approach to integrating the viscous flux is obtained by using the 
divergence theorem, 

/ V-(cc7 a F v )di1 = w 0 f"-ndT (5.86) 

Ja, Jr, 

where is the generalized control volume containing node i, with two-dimensional 
area equal to Si, and U is the boundary of IV 

Haselbacher[106] et al. have recently presented an approximate treatment for in- 
tegrating Eqn. 5.86 on two-dimensional unstructured grids, which they relate to the 
thin-layer approximation of the Navier-Stokes equations presented by Gnoffo[16] for 
structured grids. The method preserves positivity for the Laplacian and is transparent 
to grid topology. 

The development begins with the expression, 


F v -h 



f 


0 


\ 

1 


/' 

Vu-h + |V-Un x — Vv-t 



R eoo 


T 

Vv-h + \ V -VrF + Vu-t 




\nVT-h - 1 - p 

| (V-V) (V-h) + \VV 2 -h - uVv-t + vVu-t 

) 


where t is a tangent vector with components (— n y , n x ). Haselbacher’s approximation 
neglects all tangential terms and approximates V-V ~ rUVu-h + rAVv-h. Using the 
notation u n = Vu- fi. etc., and including the axisymmetric terms gives the approxi- 
mation, 


V-V —V n -h + u 7— 
V 


(5.88) 
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leading to, 

/ ° \ 

F ”-n ~ Vj + |« T (v n -h + (5.89) 

\lT n + V-V n + \VYi{Vn-n + ^)j 

A further simplification aligns n with the nearest mesh edge for faces of F on the 
interior of the domain, so that terms such as u n reduce simply to the difference in 
nodal values divided by edge length. Also, w a is chosen to be the midpoint of the 
edge. 

Including one more approximation, namely replacing the length T of the median- 
dual face with the length of the associated containment-dual face, has the effect of 
canceling some of the errors for very high-aspect-ratio cells introduced by assuming 
n is edge-aligned. For low-aspect-ratio cells, the containment dual is the same as the 
median dual and the true h is closely aligned with the mesh edge. This implementa- 
tion is similar to the suggestions of both Barth [22] and Haselbacher[107] , yet retains 
the global median-dual implementation required by a distribution scheme. 

The viscous axisymmetric source term has only one non-zero entry, 

(5 ' 90) 

Integrating using the Haselbacher approach leads to, 

I B'; d,n = I l —dM--^l ijV-ndT (5.91) 

jQi •! 'Fix, JQi V J r. ; 

Mass lumping to the node for the first term yields, 

/ — dU = Hi Si— (5.92) 

Jut y Vi 

while the second term is evaluated at edge midpoints. 

5.4 Boundary Conditions 

5.4.1 Weak Formulations 

Boundary nodes may be updated either strongly, where the nodal solution values are 
simply assigned a priori, or weakly, where the solution values at the boundary nodes 
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I = Boundary state = Reconstructed state 


Figure 5.2: Weak implementation of finite volume boundary condition for node 0, 
imposed by specifying external state. Quadrature points denoted by X ! s. 


are relaxed in time using the same formulations as for interior nodes. 

For finite volume, the weak boundary implementation specifies the solution state 
to the outside of boundary faces, then allows the approximate Riemann solver to 
construct the appropriate fluxes through the boundary faces. See Figure 5.2 for an 
illustration of the weak finite volume boundary condition. The solution state to the 
inside of the boundary face can be set from either a first- or second-order reconstruc- 
tion from the node. For some cases, second-order reconstruction to boundary faces 
has led to localized oscillations in the solution convergence at boundary nodes. 

Weak formulation of the fluctuation splitting boundary condition is developed 
using fictitious “ghost” nodes, one for each boundary node, as shown in Figure 5.3. 
Considering the scalar case, the positioning of a ghost node such that the edge con- 
necting the ghost and boundary nodes is parallel to the advection velocity results in 
a boundary fluctuation of, 

<t>\ bc 0 = — U 0 ) (5.93) 

for node 0 of Figure 5.3. Observe that this formulation is independent of the physical 
location of the ghost node, so the ghost node can be chosen to be infinitesimally close 
to the boundary node it supports. The solution state at the ghost node remains to 
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fo, fi are ficti- 
tious nodes for 
sending bound- 
ary fluctuations 
to nodes 0 resp 1. 


Figure 5.3: Weak implementation of fluctuation splitting boundary condition, im- 
posed by specifying external state at ghost nodes / 0 , f\. 


be specified, and can be varied node-to-node. The associated artificial dissipation is, 

C 0 = sign(.4-n 0 i)Co (5-94) 

and the resulting distribution is, 

S 0 U 0t t— -(Co + Co) + COE (5.95) 

Since no account of the ghost cell area is made in forming the dual area on the LHS 
of Eqn. 5.95, a scale factor on [|, 1] can be applied to the distribution. 

The extension to systems follows by analogy. The boundary fluctuation is defined, 

Co = -fCiAc 0 -noi(W /o - Wo) (5.96) 

with the artificial dissipation, 

Co = sign (Aco -foil) Co (5-97) 

Additional dissipation for the suppression of expansion shocks is added to Eqn. 5.97 
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following the form of Eqns. 5.68 and 5.69 as, 


^01 

+^a 0 — 


0 0 
0 ^oi(Aa 3 + Aa 4 ) 

0 n x 01 n y 01 ( Aas + AaJ 
0 onJ^AAs - Aa 4 ) 


0 

n 01 n 0l( A A 3 + Aa 4 ) 
n 0l(^A 3 + Aa 4 ) 
an oi ( Aa 3 - Aa 4 ) 


0 

^(Aa 3 - A a J 
^(Aa 3 - A A4 ) 
(Aa 3 + Aa 4 ) _ 


(W /o - W 0 


(5.98) 


The distribution to the boundary node is then formed as, 

•S'oUo, 4- -U(, t /(0 6co + <p bco ) + COE (5.99) 

This system treatment is only approximate, as the cross-fluctuation does not van- 
ish when V j| fo 0, as in the scalar case, but reduces to the term, 


0 0 
0 0 
0 0 
0 a 2 n% 

Strictly, there should be some cross-coupling with the neighboring boundary nodes. 
However including the term from Eqn. 5.100 requires explicitly locating the ghost 
nodes, which can be impossible for certain geometries. 

The distribution to node 1 is formed analogously, substituting 0 for 1 in Eqns. 5.93- 
5.100. 


0 

0 


“Mb 


n 


n 


(5.100) 


5.4.2 Boundary Types 

The freestream boundary condition is enacted by specifying a complete, constant 
thermodynamic state and velocity vector. By using the weak boundary enforcement, 
this one boundary condition covers the four permutations of subsonic or supersonic, 
inflow or outflow. 

The inviscid wall boundary is implemented by mirroring the primitive variables, 
either across the face for DMFDSFV or at the node for fluctuation splitting. 
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Viscous walls define a no-slip velocity and a specified wall temperature. The zero 
velocity at the wall causes the viscous axisymmetric source to be zero. 

Both the full viscous flux and the approximate thin-layer flux of Haselbacher 
reduce to, 


F v -h = 


R, 


0, pb n , nT n 


(5.101) 


at a wall, since V, V-V, and V n ■ h go to zero. Defining the heat transfer into the wall 
according to Fourier’s law, 


Qw = 


and the wall shear stress, 


t w = 

allows the Eqn. 5.101 to be written as, 


K 


R P 


-T 

T). 



v; 


(5.102) 


(5.103) 


F -77 — [0, T w . Q w ] 


(5.104) 


where the minus sign results from the choice of an outward unit normal, h, to the 
control volume, which points into the wall at a boundary. 

The solid wall is enforced weakly, by specifying the wall shear that will drive the 
flow momentum to zero and the heat flux that will drive the solution temperature to 
the desired wall temperature. An advantage of this weak approach is that wall heat 
transfer and skin friction are solved for directly, rather than as a post-processed least- 
squares reconstruction. Using an explicit update, the wall heat flux can be isolated 
as, 


Qw 


J2 r « 


w a S 

At 


(U 4 


Ui e(T w )) + RHS\ +V 


(5.105) 


Similarly, the wull shear is, 


— 




UJ n S 


At 


(u 2 , u 3 ) + rhs 1 2 % v 


(5.106) 
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On the axisymmetric axis the control- volume centroid for W a in Eqn. 5.5 is used for 
the temporal evolution, to avoid a singularity in the nodal update. The DMFDSFV 
axis condition is implemented by setting all fluxes to zero through the axisymmetric 
axis. For fluctuation splitting the axis is treated in the same manner as an inviscid 
wall, but with the radial term set to the cell centroid to ensure the numerator in 
the distribution contribution goes to zero faster than the denominator from the LHS 
time evolution. These implementations of the axisymmetric axis are observed to 
work well in practice except at the stagnation point on the axisymmetric axis, which 
is a singularity in the flowfield as well as a mathematical singularity. The deficiency 
observed with the weak boundary implementation at the axisymmetric stagnation 
point is a lack of preservation of massffow, where the velocity vectors sometimes 
allow a slight leak into or out of the domain, depending upon the particular flow 
conditions. 


5.5 Temporal Evolution 

Solutions at the nodes are updated using an explicit forward-Euler LHS. A .Jacobi 
relaxation strategy is followed with local time stepping. 

The CFL (Courant[108] et al.) criteria for explicit schemes is adapted for use with 
the node-based unstructured scheme. The inviscid timestep is defined by the most 
restrictive time for signal propagation, at the eigenvalue speeds, between adjacent 
nodes, 


I Vo •'rod + ao 1 1 rod 


where the current node is node 0 and i takes on nodal values for each distance-one 
neighbor of the current node. 

The viscous timestep restriction is taken to be an approximation based upon the 
positivity analysis for the scalar case, Eqn. 3.54, assuming order-1 Prandtl number, 


Ato 


ASoppRg^R 

Po(R + 7 — 1) J2t Jy 


(5.108) 
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The stability and convergence of axisvmmetric solutions is found to be enhanced 
by scaling the timestep for points near the cylindrical axis by the maximum of either 
the node height or the square root of the median-dual area. 

The more restrictive of the inviscid or viscous timestep is used to scale the nodal 
update. 


5.6 Verification and Validation 

5.6.1 Coding Strategy 

The code is written using a literate programming[109] style that blends source code 
(written in C) with documentation (RT£X) in the same file. Maintaining close physi- 
cal proximity between the code and documentation aids the debugging process. Mod- 
ularity is emphasized, with individual functions being verified using sample inputs 
prior to being linked with the main driver routine. 

Verification of the complete solver is performed in stages using a methodology 
derivative of Singhal[101]. A variety of canonical cases are constructed, including 
grid distortions, that are designed to exercise combinations of the various functions 
that comprise the complete solver. Validation is the application of the solver to 
complex fiowfields with comparison to benchmark data. 

5.6.2 Inviscid Verification 

Distorted mesh 

The first verification case simply passes a uniform flow through a distorted grid, with 
success being the preservation of uniformity to at least six significant digits. The 
domain is initialized to stagnant conditions with freestream flow impulsively applied 
at the boundaries. A variety of flow angles were tested on —180° < AO A < 180° for 
subsonic, transonic, and supersonic Mach numbers. Regular, high aspect ratio (100), 
skewed (2° < 9 < 175°), and randomly distorted (Figure 5.4) meshes with 121 nodes 
were used. Initial runs were instrumental in refining the treatment of eigenvalue 
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Figure 5.4: Sample randomly distorted mesh used for solver verification cases. 


limiting for fluctuation splitting as described in sections 5.2.2 and 5.4. All final runs 
were successful for both DMFDSFV and fluctuation splitting. 

Converging Mach streams 

Thermodynamic routines are verified by considering converging Mach streams, in- 
clined at ±10°. The upper stream is at Mach 2.3 while the lower stream has Mach 1.8. 
The two streams have matched densities but a temperature ratio of 1.0812, resulting 
in a horizontal slip line behind the oblique shocks. A complete description of the 
analytic solution appears in Figure 5.5 and Table 5.1. 

A sequence of four meshes, with a refinement ratio of 1.5, is considered. The 
meshes are triangulated from 16 x 16, 24 x 24, 36 x 36, and 54 x 54 grids. The 
triangulated 16 x 16 grid is shown in Figure 5.6. The finer meshes cover the same 
domain and are constructed similarly to the shown mesh. 

A Mach-number contour plot for fluctuation splitting on the finest mesh is shown 
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Figure 5.5: Description of converging Mach stream problem. Flow from left to right, 
with oblique shocks, solid, and slip-line, dashed, emanating from trailing edge of 
splitter plate. 


State 

Pi kg/nr 

! T, K 

P, kPa 

V, m/t 

A 

1.2 

300 

103.34 

798.6 

B 

1.2 

324.7 

111.73 

649.9 

C 

1.813 

356.7 

185.62 

723.7 

D 

1.718 

376.4 

185.62 

563.7 


Table 5.1: Analytic thermodynamic states for converging Mach streams. 
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Figure 5.6: 16 x 16 mesh for converging Mach streams. 


in Figure 5.7, showing crisp discontinuity resolution and the correct post-shock Mach 
numbers. The shock angles for all eight cases, i.e. DMFDSFV and fluctuation splitting 
on each mesh, are measured to be correct within ±1°. The Z, 2 -norms of the primitive- 
variables error at states C and D are plotted versus the characteristic mesh size in 
Figure 5.8. The slopes of the regression lines are indicative of the order of accuracy 
with respect to grid convergence of the two schemes for this test case. DMFDSFV 
exhibits second-order convergence, as expected. Unexpectedly, fluctuation splitting 
shows super-convergence for this particular case. True multi-dimensional upwinding 
is likely the source of the exceptional fluctuation splitting accuracy for this purely- 
supersonic flow. Supplementing the graphical determination of the grid-convergence 
rates, the equations presented by Roache[102], based on a Richardson extrapolation, 
yield average grid-convergence rates of 3.0 for fluctuation splitting and 2.1 for DM- 
FDSFV. 

Temporal convergence rates are plotted in Figure 5.9, with timings performed on 
an IRIS R10000 platform. All cases were run using the Minmod limiter and a Jacobi 
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Figure 5.8: Grid convergence rates for converging Mach stream case. Circles = fluc- 
tuation splitting, squares = DMFDSFV. 
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Figure 5.9: Convergence histories for converging Mach stream case. Fluctuation 
splitting solid, DMFDSFV clashed. Coarsest mesh on left, finest on right. 


update strategy with local time steps. Fluctuation splitting was run with a unity CFL 
number, while best convergence for DMFDSFV was found for CFL=0.7. Fluctuation 
splitting runs at 145 /is per node per iteration, while DMFDSFV" runs at 165 ps per 
node per iteration. 


Diamond airfoil 

A verification of the inviscid wall boundary condition is performed on a diamond airfoil 
at zero angle of attack and Mach 1.5. The flow deflection is five degrees. The grid 
is shown in Figure 5.10. A Mach-number contour plot using fluctuation splitting is 
shown in Figure 5.11. The corresponding DMFDSFV solution, not shown, is visually 
indistinguishable from the fluctuation splitting solution. The analytic drag coefficient, 
based on chord length, is 0.02760. The fluctuation splitting drag coefficient is 0.02638, 
for a 4.4 percent error. The DMFDSFV result has an error of 6.6 percent from a drag 
coefficient of 0.02579. 
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X 

Figure 5.10: Grid for diamond airfoil verification test. 



Figure 5.11: Mach contours on diamond airfoil, M = 1.5, fluctuation splitting solu- 
tion. 
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Figure 5.12: Two-dimensional 10 percent circular bump mesh with isobars from fluc- 
tuation splitting solution at Mach 0.1. 


Circular bump 

A subsonic two-dimensional verification is performed on a 10 percent circular bump 
at Mach 0.1. The 1389-node mesh with isobars from the fluctuation splitting solution 
is shown in Figure 5.12. A true incompressible inviscid flow would have symmetric 
isobars fore and aft, and zero drag. The fluctuation splitting drag coefficient, based 
on cord length is 0.0058. DMFDSFV predicts a drag coefficient more than twice as 
large, 0.0128. A lower fluctuation splitting drag coefficient is indicative oflower levels 
of artificial dissipation in the solution for this case. 

Sphere 

In a similar vein, Mach 0.1 axisymmetric flow over a sphere is tested on a 1369-node 
mesh. The drag coefficient, based on frontal area, is 0.43 for DMFDSFV but 0.56 for 
fluctuation splitting. Contrary to expectation, the increased artificial dissipation in 
the DMFDSFV solution creates enough of a total pressure loss to nearly eliminate 
separation on the leeside, whereas the leeside increase in pressure toward the centerline 
in the fluctuation splitting solution does produce a sizable separation region, and in 
this case a larger drag coefficient. As with subsonic bump case, true incompressible, 
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inviscid flow should theoretically produce zero drag. Both solvers converged to steady 
solutions for this case. 

Cone 

The final inviscid verification is for an 11-degree semi- vertex-angle cone at Mach 1.5. 
The well-established Taylor-Maccoll[110] method for the conical supersonic Euler 
equations predicts a drag coefficient, based on base area with no base pressure, of 
0.7795. The fluctuation splitting solution, which converged seven orders of magni- 
tude in 38 seconds, predicts a drag coefficient of 0.7785, for only a 0.13 percent error. 
The DMFDSFV solution, which took 13 percent longer at 43 seconds to reach seven 
orders of magnitude residual convergence, predicts a 0.7754 coefficient, for an error 
of 0.53 percent. 

5.6.3 Viscous Validation 

Two canonical viscous validation cases are considered: a subsonic flat plate and a 
hypersonic cylinder. Steady laminar solutions are obtained using the Haselbacher 
thin-layer viscous treatment with containment-dual modification as described in sec- 
tion 5.3. 

Flat plate 

The classic Blasiusflll] flat-plate boundary layer problem is solved on a rectangular 
domain. Mach 0.3 flow enters 2 units upstream of the plate leading edge, which is 
located at the origin. The plate is 4 units long, ending at an extrapolation outflow 
boundary. The upper boundary is 1.2 units above the plate. The Reynolds number 
is 10 4 . 

The meshes are obtained from a structured grid containing 37 equally-spaced 
points parallel to the plate, 12 points upstream of the plate and 25 points on the 
plate, and 41 points normal to the plate. The vertical grid spacing at the wall is 
0.001 units with an exponential stretching as described in Ref. [112], placing approx- 
imately 20 nodes within the boundary layer. The unstructured mesh is formed from 
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(a) Fluctuation splitting. 


(b) DMFDSFV. 


Figure 5.13: Boundary layer profiles of tangential velocity extracted from three sta- 
tions on flat plate. 


the structured grid using diagonal cuts in an alternating pattern. Two coarser meshes 
are similarly constructed by successively deleting every-other node in the wall-normal 
direction, leaving 10 and 5 nodes, respectively, in the boundary layer for the medium 
and coarse grids. 

Boundary layer profiles of u are extracted at x = 1,2,3 from both the fluctu- 
ation splitting and DMFDSFV solutions and plotted versus the Blasius solution in 
Figure 5.13. The boundary layer scaling variable is defined as, 



V = V 


(5.109) 
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(a) Containment dual approximation. (b) Strict median dual implementation. 


Figure 5.14: Boundary layer profiles computed using two different viscous dual mesh 
definitions. 


Both solution sets match the Blasius profile, indicating well-developed flow with ad- 
equate grid resolution on the finest mesh. 

Figure 5.14 shows the effect of using the containment-dual approximation in the 
Haselbacher thin-layer viscous treatment. Boundary layer profiles of u are again 
extracted at x = 1,2,3, with both solutions being run with fluctuation splitting. 
Figure 5.14(a) is the same as Figure 5.13(a), while Figure 5.14(b) uses the strict 
median-dual definition for the viscous terms. For the highly-stretched grid elements 
used in this case, it is clear that the containment-dual approximation provides im- 
proved boundary-layer resolution, while omitting the approximation leads to a profile 
that is '‘too full.” 
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Figure 5.15: Boundary layer profiles of vertical velocity extracted from midpoint of 
flat plate. 


The u-velocity profiles from the fluctuation splitting and DMFDSFV solutions 
are compared in Figure 5.15, both extracted from the plate at x = 2. The fluctu- 
ation splitting solution comes much closer to matching the Blasius profile than the 
DMFDSFV result. Excessive artificial dissipation is produced by the DMFDSFV 
scheme in the y-momentum equation, which suppresses the 'f, ’-velocity below the ana- 
lytic value. The artificial dissipation contributions to the y-momentum equation are 
plotted for both fluctuation splitting and DMFDSFV in Figure 5.16. The vertical 
scale has been enlarged by a factor of 30 to zoom in on the boundary layer in Fig- 
ure 5.16. Clearly, DMFDSFV is producing significantly more artificial dissipation 
than fluctuation splitting over the length of the boundary layer. 

For this essentially incompressible case, the suppression of the vertical velocity 
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X 

(a) Fluctuation splitting. 



- 2-10 1 2 3 4 


X 

(b) DMFDSFV. 

Figure 5.16: Artificial dissipation production in the y-momentum equation. Eleven 
contours spaced equally on 0-0.0005. 
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due to excessive artificial dissipation is manifested by an increase in skin friction co- 
efficient, as shown in Figure 5.17, where the friction coefficient increases with running 
length for DMFDSFV, but not for fluctuation splitting. Recall that DMFDSFV is 
continuously producing artificial dissipation over the length of the plate while the fluc- 
tuation splitting dissipation is restricted to the leading-edge region only. Figure 5.17 
presents data from all three grid refinement levels. The DMFDSFV results degrade 
dramatically with coarsening of the mesh, but the fluctuation splitting results remain 
relatively invariant with mesh resolution, all the way down to only five nodes in the 
boundary layer. 

The medium-mesh DMFDSFV solution was repeated using the full Navier-Stokes 
treatment, rather than the thin-layer equations. No change in the skin-friction results 
are seen over the first half of the plate, Figure 5.18, though there is an 8-percent 
improvement toward the end of the plate. Solving for the full Navier-Stokes terms 
requires 11 percent more CPU time per iteration. 

Cylinder 

The opposite end of the Mach-number spectrum is used to validate heat-transfer 
calculations, in this case for a cylinder of 1 m radius in Mach 17.6 flow. The perfect-gas 
assumption is a poor physical model for these extreme conditions, V^ = 5 km/s, Poo = 
0.001 kg/m 3 , Tqo = 200 K, T u , aM = 500 K, but the case provides a severe test of the 
algorithms under a re-entry scenario. Results are compared against the LAURA [15, 
16, 113] benchmarks. 1 The LAURA code is well-established as a structured-grid 
hypersonic solver. Also included in the LAURA benchmark data is a solution using 
the unstructured-mesh finite volume solver FUN2D[83]. The FUN2D code employs 
the same basic inviscid and viscous discretization strategy as the present DMFDSFV 
scheme. However, differences exist between the two codes in their eigenvalue limiting 
and their flux limiting, where FUN2D uses the Venkatakrishnan limiter[114] while the 
present DMFDSFV scheme is using Minmod. Also, FUN2D uses a strong boundary 
enforcement as opposed to the weak formulation and the FUN2D heating rates are 
post-processed from the temperature field rather than using the heat flux directly 

Tttp : / /hef ss . larc . nasa.gov/ 
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Figure 5.17: Skin friction coefficients for Blasius flow. 
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Figure 5.18: Effect of viscous modeling on skin friction. 


from the RHS discretization, as is used here. 

The unstructured grid for this case was obtained by simple triangulation of the 
LAURA grid, which has 65 nodes perpendicular to the surface, clustered to the wall, 
and circumferential nodes spaced every 3 degrees. Only the forward-half of the cylin- 
der is solved, as shown in the mesh and fiowfield solution of Figure 5.19. 

The surface pressure coefficient is plotted versus rotation angle from the stagna- 
tion point for both the fluctuation splitting and DMFDSFV solutions, along with the 
LAURA and FUN2D results in Figure 5.20. The LAURA, FUN2D, and fluctuation 
splitting curves all over-plot, and the DMFDSFV solution nearly over-plots, being 
1 percent low at the stagnation point and slightly high by a similar amount 90 de- 
grees away. The calculations were repeated on a grid coarsened by a factor of four 
(skip of two in both structured-grid directions), with surface pressure results plotted 
in Figure 5.21 along with the fine-mesh LAURA solution. The coarsened fluctua- 
tion splitting surface pressures retain good agreement, and the DMFDSFV solution 
matches over most of the cylinder, with minor exceptions again at the stagnation 
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Rotation degrees from stagnation point 


Figure 5.20: Cylinder surface pressures, solid = fluctuation splitting. LAURA, and 
FUN2D, while dashed = DMFDSFV. 



Rotation degrees from stagnation point 

Figure 5.21: Cylinder surface pressures on coarsened mesh, solid = LAURA, 

dashed = fluctuation splitting, and dotted = DMFDSFV. 
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Rotation degrees from stagnation point 


Figure 5.22: Cylinder surface heat-transfer rates, solid = LAURA, dashed = FUN2D, 
and dotted = fluctuation splitting. 


point, 1 percent high on this grid, and at the 90 degree point. 

Surface heat-transfer rates for LAURA, FUN2D, and fluctuation splitting are 
shown in Figure 5.22. Both of the unstructured-mesh solutions show elevated heat- 
ing at the stagnation region, with fluctuation splitting being 30 percent higher than 
LAURA while FUN2D is 50 percent higher. The DMFDSFV solution is shown in 
Figure 5.23, predicting stagnation heating rates more than double the LAURA pre- 
dictions. Also included in Figure 5.23 is the fluctuation splitting solution on the 
coarse mesh, which is seen to produce heating rates closer to the LAURA solution 
than DMFDSFV does on the fine mesh. While the basic DMFDSFV and FUN2D 
schemes are ostensibly the same on the interior domain, the significant differences in 
boundary implementation and limiting detailed at the start of this section are the 
reasons for the differing heat transfer predictions between the codes. 

The fine-mesh solutions were repeated using the full Navier-Stokes treatment, and 
no changes in heating levels were observed. 
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Rotation degrees from stagnation point 


Figure 5.23: Cylinder surface heat-transfer rates, solid = LAURA, dashed = DM- 
FDSFV(fine), and dotted = fluctuation splitting(coarse). 


5.6.4 Summary of Results 

Six inviscid test cases verify and compare the fluctuation splitting and DMFDSFV 
implementations. Both schemes are able to maintain uniform flow on a severely 
distorted mesh. DMFDSFV displays the design order of accuracy, second-order, for 
converging Mach streams while fluctuation splitting surprisingly displays third-order 
accuracy for this case. Fluctuation splitting is more accurate than DMFDSFV for 
the diamond airfoil, circular bump, and supersonic cone, but DMFDSFV is more 
accurate for the sphere test. Timings reveal fluctuation splitting runs 10-12 percent 
faster than DMFDSFV per node with the present algorithms. 

The two viscous validation cases are an incompressible flat plate and M = 17.6 
cylinder. On the flat plate fluctuation splitting produces significantly less artificial 
dissipation than DMFDSFV and provides much better skin friction predictions on 
coarse meshes, down to only five points in the boundary layer. Both schemes produce 
excellent surface pressures for the hypersonic cylinder and fluctuation splitting does 
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better for surface heating, though neither unstructured scheme predicts heating as 
well as the structured-mesh benchmark result. 



Chapter 6 

System Mesh Adaption 


6.1 Overview 

The adaption strategies detailed in chapter 4 for scalar advection-diffusion problems 
are extended to the Navier-Stokes system of equations. For the first time, the fluctu- 
ation minimization strategy for mesh adaption is developed for two-dimensional and 
axisymmetric systems, and is compared with the curvature clustering adaption that 
is representative of the current state of the art. 

This chapter’s application case is a Mach-10 wind tunnel simulation of a sting- 
mounted Mars capsule model. This case has previously been studied by Hollis[l 15] 
using a structured-grid adaption strategy, and the resulting agreement with experi- 
mental data for the sting heating was “found to be highly dependent upon grid res- 
olution and grid quality” in the wake. The structured-grid approach to this problem 
required an extremely fine mesh with extensive adaptation by hand. The challenge for 
the unstructured-grid techniques will be to obtain comparable accuracy using fewer 
grid points with fully automated adaption. 


6.2 Curvature Clustering 

The local, anisotropic adaption strategy based upon a posteriori error estimates is 
extended directly from section 4.1. In extending to the system, Eqn. 4.2 generalizes 
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A scalar value for the error estimate is required, with a simple choice being to just 
take the Z^-norm of the error vector. Also, the error does not necessarily need to 
be formed from the conserved variables, but could alternatively be formed from the 
primitive variables, the Mach number as Haba.shi[76] does, or even a derived quantity 
such as heat transfer. 

The basic operations remain point deletion, edge swapping, point insertion, and 
nodal displacements. These procedures are applied sequentially in the order listed, 
with the solver iterating on the solution between each of the four steps, constituting 
a complete adaption cycle. A tolerance can be set for thresholds to define when the 
mesh is sufficiently adapted. 

While the explicit use of nodal gradients in Eqn. 6.1 makes this style of adaption 
naturally symbiotic with the DMFDSFV scheme used here, there is no restriction 
against using curvature clustering with the fluctuation splitting solver. Also note 
that Eqn. 6.1 is written in vector notation, and so can be applied directly in three 
dimensions just as it is used here in two spatial dimensions. 


6.3 Characteristic Alignment 

Chapter 4 demonstrated that fluctuation splitting is an exact solver for linear ad- 
vection on a characteristically-aligned mesh. Adaption for characteristic alignment 
was developed as an automated, anisotropic local-operation strategy in sections 4.2 
and 4.4. Applications to non-linear and advection-diffusion problems showed that 
alignment with characteristics was approximated, but not fully achieved, by the au- 
tomated local scheme. For the Navier-Stokes system, true characteristic alignment is 
not a physically valid ideal due to the presence of multiple and/or imaginary charac- 
teristics. However, the automated local anisotropic adaption strategy is extended here 
with the objective of achieving solution improvement with fluctuation minimization 
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as the agent to obtain characteristic-driven alignment in supersonic regions consistent 
with mesh enrichment in the subsonic regions. 

As with the scalar case, the analogy to the edge error estimate is formed from 
the fluctuations, both inviscid and viscous, in the cells adjoining the edge. For ax- 
isymmetric cases the fluctuation is scaled by the inverse of the cell-centroid y value. 
Additional weighting by the inverse of the square root of the cell area is found to 
balance the contributions from neighboring cells of disparate sizes. 


Point deletion 

Nodes are flagged for deletion if the sum of the L 2 -norms of the fluctuations in all 
surrounding cells is below a threshold. Edges are swapped in an attempt to reduce 
the number or edges connected at the node to three or four to match the canonical 
point-removal patterns. If the local connectivity is too complicated for the automated 
pattern matching, the node is simply left in the domain. During a typical cycle, on 
the order of only 10 percent of the nodes that had been flagged for deletion will be 
left in the domain because of localized complicated connectivity. 


Point insertion 

An edge is split, adding a node, if the sum of the L 2 -norms of the fluctuations in the 
cells to either side of the edge exceed a threshold. The cell fluctuations are weighted 
in the same manner as was described for the point deletions. 


Edge swapping 

An edge is flagged as a candidate for swapping if the RMS of the fluctuations in the 
cells to either side of the edge exceed the target threshold. If the swapped config- 
uration maintains a physically valid grid, then the RMS of the fluctuations in the 
swapped cells are computed. If the swapped configuration has a smaller RMS value 
then the edge swap is performed. 
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Nodal displacement 

Mesh movement for the system is driven by the minimization of a functional, just as 
for the scalar case. The system extension of Eqn. 4.21 is, 

Tj =2 (6-2) 

T 


The functional T can be expressed in terms of either the conserved or auxiliary 
fluctuations, as the minimization of one implies the minimization of the other. The 
derivative of T (see Eqn. 4.22) is formed using the chain rule as, 


ar* 

dxi 





1 , H ^ s t 


0 T + 0 T S T 



(6.3) 


Having defined the gradient of the functional, the method of steepest descent from 
Eqns. 4.24 and 4.25 can be applied directly to drive the nodal displacements. 

The weighting factor Sj is a symmetric positive-definite matrix. Weighting each 
equation equally without regard to cell sizes results in Sj = /, while inverse area 
weighting yields Sj = jrl. The derivatives of S T for these and other area-weighted 
choices have been covered in Eqns. 4.50-4.56. 

The cell fluctuation can be rearranged from Eqn. 5.41 (two-dimensional terms) as, 


0 


3 = 1 


(6.4) 


where W j = W^Z,. Following Eqns. 4.32 and 4.33, the derivatives of the fluctuation 
are, 
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The flux Jacobian of the auxiliary variables is. 
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which leads to the approximation ~ | yl. As was done for the scalar case, it is 
assumed that moving a node does not change the solution at the other two nodes of 
the triangle, so that the variation of the cell-average flux Jacobian scales like one-third 
the variation of the velocity at the node being moved, fy ~ 

Similarly, the variation of the Jacobian of the transformation scales like one-third 
the nodal variation, ~ and is neglected as a sub-principle term relative to 
the change in solution value, 


awj 

dx 2 


W 2 ^ 

dx 2 


( 6 . 8 ) 


Since the solution is locally assumed to vary only at the node being moved, i.e. , 

dZ i _ ^ _ g 
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3 =1 


(6.9) 


The remaining term to evaluate is yy. The steady-state distribution can be 
written, 


53 ^ + M a ) a - (7 + M^) p\ W 2 Z 2 


= ^[(/ + M tt )aW^ 
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(/ + MtiPYrM ( 6 . 10 ) 


or, 
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( 6 . 11 ) 
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where, 

a + = — (/ + M a ) a (6.12) 

/ 3 + =\(I + M P )(3 (6.13) 

Differentiating while freezing the Jacobians leads to, 






The derivatives of a and /3 remain as Eqns. 4.41 and 4.42. 

Having determined ^ follows from, 

d = d{y/pu) -ud(y/p) dv = d{y/pv) -vd{y/p) 19 

u \fp v \fp 

The concept of fluctuation minimization extends directly to three dimensions for 
point deletion and insertion. The idea of edge swapping extends as face reconnec- 
tions amongst tetrahedra in three dimensions which is more complicated due to the 
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increased degrees of freedom and geometric connection possibilities. The functional 
for nodal displacements. Eqn. 6.2, is written in vector notation and so can be applied 
directly to three dimensions. However, not only would the derivatives of the fluc- 
tuation splitting scheme need to be reformulated, but the basic fluctuation splitting 
scheme itself would need to be worked out in detail to move to three dimensions. 

A three-dimensional alignment strategy has been employed by Beeman and Pow- 
ers[l 16] . That method starts on an attached shock and projects downstream, trying to 
remain on the shock. This involves iterating an assumed shock shape and an implied 
body, iterating until the implied body matches the true surface. Open questions 
remain on how to start on detached shocks and how to handle embedded shocks. A 
significant drawback mentioned by the authors is that, “[The method] far enough 
downstream will eventually become unstable.” 


6.4 Mars Pathfinder 

The demonstration case for the system adaption is borrowed from the Mach-10 wind 
tunnel tests of Hollis[l 15, 117, 118], which investigated the aerothermodynamic envi- 
ronment experienced by a payload in the wake of an aerobrake. For his dissertation, 
Hollis specifically looked at the effect of grid adaption on sting heating for his model. 
Using a structured-grid finite volume solver, a mesh of 125 x 357 = 44,625 nodes was 
required to get reasonable pressure and heating-rate agreement on the sting. Addi- 
tionally, this fine mesh required extensive and time consuming adaption by hand to 
sufficiently resolve the wake flow. 

Using the Hollis data as a benchmark, the present study seeks to achieve compara- 
ble results on an automatically-adapted unstructured mesh. Significant time savings 
can be achieved primarily by automating the adaption procedure, but also reducing 
the solver iteration and convergence times if the unstructured approach requires fewer 
grid nodes for the same resolution. 
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6.4.1 Configuration and Conditions 

The axisymmetric wind tunnel model of the Mars Pathfinder capsule consists of a 
spherically-blunted 70-degree sphere-cone. The radius is 1 in. and the nose radius is 
\ in. The shoulder radius is ^ in. and the aftbody angle is 40 degrees. The base 
radius is | in. The sting is 4 in. long with a radius of || in. Figure 6.1 shows the 
boundary of the computational domain. 

The freestream conditions for the NASA Langley 31-Inch Mach 10 Air Tunnel, 
corresponding to a nominal Reynolds number per foot of 0.5 x 10 6 ft -1 , are: P = 69 Pa, 
T = 53 K, p = 0.0045 kg/m 3 , and u = 1416 m/s. The wall temperature is taken 
to be a uniform 300 K. Historical experience with this tunnel at these conditions 
indicate that laminar perfect-gas calculations are adequate for comparison with the 
experimental data. 
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In addition to the axisymmetric configuration, the present study also considers 
a two-dimensional version of the problem, using the same geometry and freestream 
conditions. 


6.4.2 Benchmark Data 

For the axisymmetric configuration Hollis has provided the experimental heat-transfer 
data along with numerical heat-transfer and surface pressure results. The numerical 
results were generated using the NEQ2D code of Candler[119]. A second set of numer- 
ical results were obtained for the present study using the LAURA code of Gnoffo[15]. 
Both of these codes are structured-grid finite volume solvers. Surface pressure co- 
efficients from the numerical results are shown in Figure 6.2 and heat-transfers for 
all three datasets are in Figure 6.3, where the size of the symbols is indicative of 
the uncertainty in the experimental data reported by Hollis, ±7%. The forebody 
pressures are in agreement between NEQ2D and LAURA, but the aftbody pressure 
agreement is weaker. The computed heating rates match the experimental data, aside 
from an irregularity in the LAURA result at the stagnation point, along the body to 
a running length of 3 inches, which is located on the sting. Hollis speculates that the 
free shear layer in the wake may be transitioning at the point where the experimental 
and computational heating rates diverge on the sting. 

For the two-dimensional configuration, the axisymmetric grid was enlarged to be 
125 x 513 = 64, 125 nodes so as to capture the bow shock due to the greater shock 
stand-off distance. A LAURA solution was obtained for this case, and the results of 
grid convergence in the body-normal direction on pressure and heating are shown in 
Figures 6.4 and 6.5. Surface pressures are largely converged with 129 points in the 
body-normal direction, and certainly by 257 points. A similar statement applies for 
the heating, except in the vicinity of the stagnation point. The heating prediction at 
the nose continually increases with each grid doubling. 
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Figure 

dashed 



(a) Heatshield. 



(b) Aftbody and sting. 


6.2: Axisymmetric capsule benchmark surface pressures, solid=NEQ2D, 

=LAURA. * 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.3: Axisymmetric capsule benchmark surface heating, symbols=experiment 
solid=NEQ2D, dashed=LAURA. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.4: Two-dimensional capsule benchmark surface pressures, solid=513 

dashed=257, dotted=129, dash-dot=65 points normal to the surface. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.5: Two-dimensional capsule benchmark surface heating, solid=513, 

dashed=257, dotted=129, dash-dot=65 points normal to the surface. 
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6.4.3 Unadapted Baseline 

Prior to performing any adaption, both unstructured DMFDSFV and fluctuation 
splitting schemes were used to obtain solutions on triangulated versions of the struc- 
tured meshes, as an indication of a target accuracy for the adaption strategies. During 
this process, a number of problems and sensitivities were encountered. 

As a general note the explicit single-stage time integration scheme used here is 
severely restricted in the stable timestep for high-Reynolds-number viscous problems. 
A structured-mesh solver, such as the LAURA code, would typically alleviate this 
stiffness by employing an implicit solver. Another option would be to use a parallelized 
strategy. However, both the development of an unstructured implicit solver for the 
present schemes and the parallelization of the schemes were considered beyond the 
scope of the current work. As a consequence, viscous solutions for the capsule on the 
fine meshes required weeks of CPU time on a desktop PC. 

Since both the DMFDSFV and the fluctuation splitting schemes are implemented 
at the nodes, the timestep on the axisymmetric axis vanishes, slowing temporal con- 
vergence. The radial term in the timestep determination is chosen at the median-dual 
centroid for nodes on the axisymmetric axis, rather than at y = 0, allowing for a more 
realistic wave propagation on the axis. The solutions are observed to remain tempo- 
rally stable with this implementation while yielding much faster convergence rates. 
Yet, the axisymmetric axis is always observed to converge slower than the rest of the 
fiowfield, due to the radial weighting on the control volume. In general, axisymmetric 
solutions take two to three times longer to converge than an equivalent two-dimen- 
sional case because of this axis singularity effect. 

The carbuncle effect, Figure 6.6, reared its ugly head for some cases, typically on 
meshes that are fine in the streamwise direction but very coarse in the body-normal 
direction. Although the carbuncle effect is a known deficiency 1 in the structured 

1 The carbuncle effect was first reported by Peery and Imlay[120] for a Mach-6 cylinder as an 
unexplainable “protuberance” in the bow shock at the symmetry plane. Later work by Roberts[121] 
and Quirk[122] reveal the mathematical basis for the production of spurious solutions for grid-aligned 
shocks with the Roe flux difference splitting scheme. As Quirk reports, “[Pjarallel to the shock, Roe’s 
scheme will not add any dissipation via the contact and shear waves, to counteract perturbations 
that appear through the acoustic waves; this appears to be a recurring theme whenever Roe’s method 
fails. R is interesting to note that if Harten’s entropy fix is applied to the contact and shear waves, 
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Figure 6.6: Fluctuation splitting u-velocity contours, showing carbuncle formation 
for two-dimensional capsule case. 

Roe solver, its behavior in unstructured schemes has not been as well characterized. 
Lowering the CFL number by three orders of magnitude did not change the solution, 
verifying that the phenomenon was not instability induced. Similarly, re-running 
time-accurately produced the same behavior, strongly suggesting the phenomenon 
was not transient induced. The traditional fix for the Roe scheme involves adding 
dissipation through eigenvalue limiting. The present schemes were being run with 
thin-layer viscous terms and eigenvalue limiting on V ± a only. Switching to full 
viscous terms and limiting all eigenvalues both served to reduce the occurrence of the 
carbuncle effect, Figure 6.7. 

The last problem encountered was that the surface heating at the stagnation point 
had the wrong trend, dropping toward zero heating. Figure 6.8 shows the f'orebody 

any shortcoming of Roe’s scheme is invariably cured. However, there is no justification, either 
physical or mathematical, for applying this fix to these waves, it is just a convenient method for 
introducing an amount of artificial dissipation into the scheme.” 
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M oo = l0 



Figure 6.7: Fluctuation splitting u-velocity contours for two-dimensional capsule case. 


heating on a two-dimensional 63 x 129 mesh using fluctuation splitting. Note the 
dramatic drop in heating at the stagnation point, contrasted with the benchmark 
data of Figure 6.5(a). Looking into the fiowfield for an explanation of the heating 
trend a recirculation bubble is found at the nose, Figure 6.9, which physically should 
not be present. The cause of this stagnation recirculation proved extremely time 
consuming to identify and it was unexpectedly found to be a grid-induced feature. 
The unstructured meshes for the unadapted cases were initially created by a simple 
triangulation of the structured grids. All quadrilateral cells were cut in the same 
direction, and this bias in the diagonals along the axis boundary produced a shift 
in the stagnation point off the axis. Repeating the calculations on a biased grid 
covering both the upper and lower half-planes removed the vortex but still shifted 
the stagnation point. The cure for this problem was to remove the grid bias by 
alternating the diagonals in the derived unstructured meshes, yielding the vortex-free 
solution in Figure 6.10 for the same case as is shown in Figure 6.9. 
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Figure 6.8: Incorrect forebody heating trend, two-dimensional fluctuation splitting. 


Y, in. 



X, in. 


Figure 6.9: Stagnation point vortex. 
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Figure 6.10: Correct streamlines in stagnation region. 


32 x 65 

32 x 129 


63 x 65 

63 x 129 

63 x 257 

125 x 65 

125 x 129 

125 x 257 125 x 513 


Table 6.1: Matrix of grid dimensions for two-dimensional baseline cases. 

Two-dimensional 

Viscous solutions for the two-dimensional capsule were computed using both fluctu- 
ation splitting and DMFDSFV discretizations on a sequence of unadapted grids as 
listed in Table 6.1. 

Results using only 32 points to define the surface were generally under-resolved for 
both schemes. Results for 63 and 125 surface points are compared in Figures 6.11-6.14 
using 257 points normal to the body. 

DMFDSFV surface pressures, Figure 6.11, overlap the benchmark results on the 
forebodv for both grids. Aftbody pressure agreement is poor for 63 surface points. 
The solution improves by going to 125 points, though the pressure coefficient is still 
over-predicted by a factor of 2. The heating trends for DMFDSFV are reversed 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.11: Two-dimensional capsule surface pressures: ^-refinement with DMFDS 
FV, solid=63, dashed=125, dash-dot-dot =benchmark. 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.12: Two-dimensional capsule surface heating: z-refinement with DMFDSFV, 
solid=63, dashed=125, dash-dot-dot=benchnrark. 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 

Figure 6.13: Two-dimensional capsule surface pressures: ^-refinement with fluctuation 
splitting, solid=63, dashed=125, dash-dot-dot=benchmark. 
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(a) Heatshield. 





_l 1 1 1 I I I I I I I I I I I I I I 

2 3 4 5 


s, in. 


(b) All, body and sting. 


Figure 6.14: Two-dimensional capsule surface heating: i-refinement with fluctuation 
splitting, solid=63, dashed=125, dash-dot-dot=benchmark. 
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from the pressure trends. Figure 6.12 shows that the coarser-mesh forebody heating 
is in poor agreement with the benchmark. The finer mesh improves the heating 
prediction, but is still more than 20 percent higher than the benchmark over much 
of the heatshield, especially at the shoulder. The aftbody and sting heating rates are 
reasonable on both meshes. 

Fluctuation splitting surface pressures, Figure 6.13, are in very good agreement 
with the benchmark data, both fore and aft. The 125 surface point mesh is in par- 
ticularly good agreement, over-plotting the benchmark over all but the tail end of 
the sting. Fluctuation splitting surface heating, Figure 6.14, is in good agreement 
for both meshes, with the finer mesh providing better resolution at the shoulder. 
The glaring weakness in the fluctuation splitting solutions is at the stagnation point, 
where the heating spikes 40 percent higher than the benchmark. Further, the known 
correct heating trend at a stagnation point calls for the heating to level off, as the 
benchmark results do, rather than spike, as with the present results. 

Grid convergence trends in the body-normal direction, for 125 surface points, are 
shown in Figures 6.15-6.18. 

The DMFDSFV surface pressure grid convergence is shown in Figure 6.15. On 
the forebody excellent agreement with the benchmark is seen for all but the coarsest 
mesh. On the aftbody and sting, only the solution from the mesh with 129 points 
normal to the body is in good agreement with the benchmark. The trend with 
refining the mesh is toward increasing pressures on the aftbody, and the solutions do 
not appear to be converging toward the benchmark result. The DMFDSFV heating 
with grid refinement is shown in Figure 6.16. On the heatshield the solution from the 
j=129 mesh provides the closest match to the benchmark, while the 257-point result 
is not too much worse. The solution from the coarsest mesh is wildly unresolved and 
the result from the finest mesh is disappointing. On the aftbody, again the coarsest 
mesh is under- resolved. The three finer meshes do appear to be converging with grid 
refinement and have reasonable agreement with the benchmark. 

The fluctuation splitting surface pressures, Figure 6.17, show good grid conver- 
gence trends, with both of the two finest meshes matching the benchmark over the 
entire body, and none of the grids produced terrible results. Fluctuation splitting 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.15: Two-dimensional capsule surface pressures: j-refinement with DMFDS- 
FV, solid=513, dashed=257, dotted=129, dash-dot=65, dash-dot-dot=benchmark. 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.16: Two-dimensional capsule surface heating: j-refinement with DMFDSFV, 
solid =51 3, dashed=257, dotted=129, dash-dot=65, dash-dot-dot=benchmark. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.17: Two-dimensional capsule surface pressures: j-refinement with fluc- 
tuation splitting, solid=513, dashed=257, dotted=129, dash-dot=65, dash-dot 
dot-benchmark. 
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32 x 45 

32 x 89 


63 x 45 

63 x 89 

63 x 177 

125 x 45 

125 x 89 

125 x 177 125 x 353 


Table 6.2: Matrix of grid dimensions for axisyrnmetric baseline cases. 

heatshield heating, Figure 6.18, also shows good grid convergence trends, with both 
of the two finest-mesh solutions matching the benchmark well except at the stagna- 
tion point, where the 513-point results do better but still have a spike of 15 percent 
at the symmetry line. Heating results in the wake region also show some grid conver- 
gence trends, though the finest mesh resolves an additional wake vortex on the sting, 
producing an inflection in the data. The coarser meshes actually match the bench- 
mark better on the aftbody, before elevating on the sting where the finer meshes, and 
the 257-point result in particular, do a better job of matching the benchmark. 

The two-dimensional unadapted baseline results are chosen to be the solutions on 
the triangulated 125 x 257 mesh. The fluctuation splitting solutions on this mesh 
were as good as on the finest mesh but with half the points. For DMFDSFV the 
125 x 129 mesh produced results generally as good as or better than the baseline, 
though the results are difficult to interpret due to the lack of clear grid-convergence 
trends. The baseline surface pressures and heating rates are shown in Figures 6.19 
and 6.20. Both solvers match forebody pressure with the benchmark. Fluctuation 
splitting also matches the benchmark pressure on the sting, while the DMFDSFV 
baseline over-predicts the pressure in this region. The fluctuation splitting forebody 
heating matches the benchmark well except as previously discussed at the stagnation 
point. The DMFDSFV forebody agreement is not as good. However, both schemes 
yield comparable aftbody heating rates. 

Axisyrnmetric 

Viscous solutions for the axisyrnmetric capsule configuration were obtained on a se- 
quence of unadapted meshes formed as triangulations from the structured meshes 
used for the benchmark computations. The sequence of grids is listed in Table 6.2. 
All of the axisyrnmetric solutions to be shown suffer from irregular heat transfer 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.18: Two-dimensional capsule surface heating: j-refinement with fluc- 

tuation splitting, solid=513, dashed=257, dotted=129, dash-dot=65, dash-dot- 
dot=benchmark. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.19: Two-dimensional baseline capsule surface pressures, solid=fluctuation 
splitting, dashed=DMFDSFV, dash-dot-dot=benclinrark. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.20: Two-dimensional baseline capsule surface heating, solid=fluctuation 
splitting, dashed=DMFDSFV, dash-dot-dot=benclinrark. 
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patterns at the axis stagnation point, true for both the DMFDSFV and fluctuation 
splitting formulations as here implemented. The stagnation point surface pressures, 
however, are well behaved. The difficulty with the heating appears to be a result of 
applying the weak form of the boundary conditions at both a mathematical singular- 
ity, the axisymmetric axis, and a surface flow singularity, the stagnation point, for the 
node-based schemes. Cell-based, rather than node-based, schemes or strong bound- 
ary enforcement formulations might not experience similar stagnation point heating 
difficulties for the axisymmetric cases. These heating irregularities are confined to 
the stagnation node and immediate neighbor nodes. 

As with the two-dimensional results, the axisymmetric solutions using only 32 
points to define the surface were under resolved. Results using 63 and 125 points 
to define the surface are presented in Figures 6.21-6.24, all with 89 points in the 
body-normal direction. 

DMFDSFV surface pressures, Figure 6.21, overlap the benchmark results on the 
f'orebody, especially on the finer mesh. Aftbody pressures are poorly resolved on the 
coarser mesh. The finer mesh reasonably matches the benchmark except right at the 
shoulder where the overexpansion is missed. DMFDSFV heating trends with stream- 
wise refinement are shown in Figure 6.22. The forebody heating is much improved 
on the finer mesh, though still 20 percent higher than the benchmark over most of 
the heatshield. The stagnation point heating is over-predicted by as much as 40 per- 
cent on both meshes and the peak heating on the shoulder is greatly over-predicted. 
However, on the aftbody and sting the finer mesh results are in very good agreement 
with the benchmark results. The coarser mesh under-predicts the wake heating. 

Fluctuation splitting surface pressures for streamwise grid refinement are shown 
in Figure 6.23. The forebody pressures are the same for both meshes, but are 3-5 per- 
cent lower than the benchmark, except at the stagnation point where the benchmark 
is matched. Aftbody pressures match the benchmark extremely well for both grids. 
Forebody heating with fluctuation splitting, Figure 6.24, matches the benchmark 
fairly well except for over-predicting the shoulder heating spike and the stagnation 
point maximum. The finer mesh does a better job at the shoulder but results from 
neither mesh are particularly encouraging at the stagnation point. The fluctuation 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.21: Axisymmetric capsule surface pressures: ^-refinement with DMFDSFV, 
solid=63, dashed=125, long-dash=NEQ2D, dash-dot-dot=LAURA. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.22: Axisymmetric capsule surface heating: ^-refinement with DM- 

FDSFV, solid=63, dashed=125, long-dash=NEQ2D, dash-dot-dot=LAURA, cir- 
cles=experiment. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.23: Axisymmetric capsule surface pressures: ^-refinement with fluctuation 
splitting, solid=63, dashed=125, long-clash=NEQ2D, dash-dot-dot=LAURA. 
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splitting heating results on both meshes provide excellent matches with the bench- 
mark data on the aftbody and sting. 

Grid convergence trends in the body-normal direction, for 125 surface points, are 
shown in Figures 6.25-6.28. 

DMFDSFV forebody pressures match the benchmark on all grids, Figure 6.25. 
Flowever, the aftbody and sting pressures do not show grid convergence. The best 
match with the benchmark comes on the grid with 89 points normal to the body. It 
is disappointing that the pressures deteriorate on the finer meshes, with the details 
being washed out perhaps by a different separation pattern in the wake. The DM- 
FDSFV forebody heating, Figure 6.26, also fails to show a good grid convergence 
trend, with the results deteriorating on the finer meshes. The two coarser meshes 
provide the best agreement, with the 89-point results being smoother than the 45- 
point solution. Still, the forebody heating is over-predicted by 20 percent over most 
of the heatshield and by more at the stagnation point. On the aftbody and sting only 
the 89-point solution provides good agreement with the benchmark data. 

Surface pressure trends with grid refinement in the body-normal direction using 
fluctuation splitting are shown in Figure 6.27. The forebody pressure is grid con- 
verged on all meshes, but slightly under-predicts the benchmark results. Aftbody 
and sting pressures on all meshes show excellent agreement with the benchmark solu- 
tions. Forebody heating with fluctuation splitting, Figure 6.28, shows good agreement 
over most of the heatshield for all but the 45-point mesh, which itself is only 20 per- 
cent elevated. As usual the stagnation point heating is largely in error. Heating in 
the wake shows good agreement for all but the finest mesh, where the heating on the 
sting shows an unexpected inflection point. 

The axisymmetric baseline solutions are chosen to be the results on the triangu- 
lated 125 x 89 mesh. The fluctuation splitting solutions are consistently grid con- 
verged on this mesh, and the DMFDSFV results on the mesh are generally the best, 
often better even than on the finer meshes. The axisymmetric baseline pressures are 
presented in Figure 6.29. On the heatshield the DMFDSFV solution overlaps the 
benchmark pressures, while the fluctuation splitting result is steadily 3 percent low. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.24: Axisymmetric capsule surface heating: ^-refinement with fluctuation 
splitting, solid=63, dashed=125, long-dash=NEQ2D, dash-dot-dot=LAURA, cir- 
cles=experiment. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.25: Axisymmetric capsule surface pressures: j-refinement with DMFDS- 
FV, solid=353, dashed=177, dotted=89, dash-dot=45, long-dash=NEQ2D, dash-dot- 
dot—LAURA. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.26: Axisvmmetric capsule surface heating: y -refinement with DMFDS- 
FV, solicl=353, dashed=177, dotted=89, dash-dot=45, long-dash=NEQ2D, dash-dot- 
dot=LAURA, circles=experiment. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.27: Axisymmetric capsule surface pressures: j-reflnement with fluctua- 
tion splitting, solid=353, dashed=177, dotted=89, dash-dot=45, long-dash=NEQ2D, 
dash-dot-dot=LAURA. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.28: Axisymmetric capsule surface heating: ^-refinement with fluctua- 

tion splitting, solid=353, dashed=177, dotted=89, dash-dot=45, long-dash=NEQ2D, 
dash-dot-dot=LAURA, circles=experiment. 
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On the af'tbody and sting it is the fluctuation splitting solution that provides an ex- 
cellent match to the benchmark. The DMFDSFV solution is generally good except 
for missing the overexpansion pressure drop at the shoulder, seen as the pressure min- 
imum at s = 1.2. The baseline surface heat transfer rates are shown in Figure 6.30. 
The fluctuation splitting heating on the forebody is in very good agreement with 
the benchmark data except at the stagnation point, where the match is poor. The 
DMFDSFV result generally over-predicts the forebody heating, except for a region 
between s = 0. 2-0.4 near the nose, where the agreement with the benchmark data 
is very good. Heating levels at the stagnation point and shoulder are both greatly 
over-predicted. Both schemes match the benchmark heating datasets on the af’tbody 
and sting. 

6.4.4 Solution- Adapted Results 

The adaption mechanisms, point deletion, edge swapping, nodal displacements, and 
point insertion, are applied one at a time, using both the curvature clustering with 
DMFDSFV and the fluctuation minimization with fluctuation splitting, to assess the 
effect of each individual component. Then a full adaption cycle is applied using a 
strategy based upon the results of the individual adaption tests. 

Two-dimensional 

For point deletion the starting solution is the converged baseline. The baseline mesh 
has 32,125 nodes, and the objective is to remove points while maintaining the accuracy 
of the baseline solution. Both schemes were able to remove 10 percent of the nodes 
(3284 nodes for fluctuation minimization and 3366 for curvature clustering) causing 
minimal change to the solution. The nodes were predominantly removed from the 
freestream. A further 10 percent of the nodes were removed (3157 for fluctuation min- 
imization and 3076 for curvature clustering), for a total of 20 percent of the nodes 
deleted from the initial grid. The heating rates for both schemes are essentially un- 
changed, Figures 6.31 and 6.32, but both schemes have trouble maintaining the bow 
shock capture due to the further loss of freestream points, as shown in Figure 6.33 for 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.29: Axisymmetric baseline capsule surface pressures, solid=fluctuation split 
ting, dashed=DMFDSFV, long-dash=NEQ2D, dash-dot-dot=LAURA. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.30: Axisymmetric baseline capsule surface heating, solid=fiuctuation 

splitting, dashed=DMFDSFV, long-dash=NEQ2D, dash-dot-dot=LAURA, cir- 
cles=experiment. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.31: Two-dimensional capsule surface heating after coarsening with curvature 
clustering, solid=baseline, dashed=10% removed, dotted=20% removed. 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.32: Two-dimensional capsule surface heating after coarsening with fiuctua 
tion minimization, solid=baseline, dashed=10% removed, dotted=20% removed. 
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Figure 6.33: Loss of shock capture due to coarsening, fluctuation splitting. 


fluctuation splitting with fluctuation minimization. A similar problem also occurred 
using DMFDSFV with curvature clustering. The removal of more points causes rapid 
solution deterioration due to the failure of both schemes to property capture the bow 
shock on the coarsened meshes. 

For evaluating the edge swapping adaption the starting mesh has half as many 
points as the baseline, about 16,000 nodes and 50,000 edges, with the goal of realizing 
an improvement in the solution. Curvature clustering was used to swap the 538 
edges (approximately 1 percent) with the largest error estimate. These edges were all 
located in vicinity of the bow shock. The resulting heat transfer rates show only minor 
changes from the starting solution, Figure 6.34, without showing a clear progression 
toward the benchmark dataset. Edge swapping with fluctuation minimization is not 
beneficial for this case. About 1 percent of the edges were swapped but the resulting 
fluctuation splitting solution does not converge due to ringing of the bow shock near 
the stagnation point. Figure 6.35 shows pressure contours in the stagnation point 
region. Notice the irregular contours downstream of the shock. The difficulty in this 
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(a) Heatshield. 
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(b) Aftbody and sting. 


Figure 6.34: Two-dimensional capsule surface heating after swapping with curvature 
clustering, solid=start, dashed=l% swapped, dash-dot-dot=benchmark. 
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Figure 6.35: Fluctuation splitting shock bulging at stagnation point after edge swaps, 
pressure contours. 

solution is that the bow shock is making discrete jumps across grid lines, induced by 
swapped edges at the shock. The unadapted bow shock more closely follows the grid 
lines at the y = 0 symmetry plane. Recall that this mesh was originally well-aligned 
with the bow shock from the structured-grid solvers, and the effect of edge swapping 
on a more random unstructured mesh may be different. Also, on a finer mesh, such as 
is used for the baseline solutions, the smaller grid spacing at the shock might lessen 
the detrimental effects of discrete jumps in the shock location. Edge swapping was 
tried again for fluctuation minimization, this time only swapping half a percent of 
the edges, but still produces the same disappointing results. 

Nodal displacements are tested using the same 16,000-node mesh as was used to 
start the edge swapping, also with the goal of improving the solution toward the 
baseline. Curvature clustering w r as used to move 1717 nodes with the largest error 
estimates, about 10 percent of the nodes. The total distance moved was 10.34 in., the 
average distance moved was 0.006 in., and the RMS distance moved was 0.0095 in. 
Nodes were moved in the shock, at the shoulder, and in the forebody boundary layer. 
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Figure 6.36: Nodal displacements at shock using curvature clustering, original mesh 
shaded. 

A close-up of movement at the bow shock is shown in Figure 6.36, where the original 
mesh is shaded. Observe that the movement is not aggressive, but rather produces 
a gentle clustering toward the shock. Surface heating rates for this case are shown 
in Figure 6.37, where a minimal change in the solution is seen. Displacing nodes 
with fluctuation minimization moved 1827 nodes a total distance of 0.7488 in., for an 
average of 0.0004 in. and RMS of 0.0015 in. The fluctuation minimization movement, 
while moving a few more nodes than the curvature clustering, moved the nodes much 
smaller distances. Again, the movement occurred at the shock, near the shoulder, 
and in the forebody boundary layer. Surface heating rates for this case are shown 
in Figure 6.38. The solution has been made much worse on the forebodj^, while not 
much change is seen in the wake. Investigating the solution reveals that a windside 
vortex pattern has emerged, Figure 6.39, that is causing the lower heat transfer rates 
on the heatshield. The solution was run a further 800,000 iterations with no change 
to this windside vortex pattern. 

Point insertion is tested on the unadapted 16,125-node mesh, seeking to add 
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(a) Heatshield. 
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(b) Aftbody and sting. 


Figure 6.37: Two-dimensional capsule surface heating after moving with curvature 
clustering, solid=start, dashed=10% moved, dash-dot-dot=benchmark. 
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(a) Heatshield. 
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(b) Aftbody and sting. 


Figure 6.38: Two-dimensional capsule surface heating after moving with fluctuation 
minimization, solid=start, dashed=10% moved, dash-dot-dot=benchmark. 
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Figure 6.39: Windside vortices produced by fluctuation minimization nodal displace- 
ments. 

20 percent more nodes. Curvature clustering added 3154 nodes, producing the heating 
results shown in Figure 6.40. On the forebody, there is improvement in the stagnation 
point region, but worsening near the shoulder. On the aftbody and sting the heating 
moves somewhat away from the benchmark result. Fluctuation minimization added 
3299 nodes, with heating results shown in Figure 6.41. On the heatshield the heating 
levels rise markedly toward the benchmark result, but exhibit a high-frequency oscil- 
lation of significant amplitude. There is not much change in the heating in the wake 
region. 

For the full adaption the starting solution is taken on the triangulated 125 x 129 
mesh with 16,125 nodes, half the number of nodes used to generate the baseline 
unadapted solutions. The intention is to look for an improvement in the coarse mesh 
solution without increasing the number of nodes. Based upon the results of the 
component adaption tests, the strategy for an adaption cycle is to delete 10 percent 
of the nodes, swap 1 percent of the edges, move 5 percent of the nodes, and then 
insert back in 10 percent of the nodes. The solution is re-converged between each 
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(a) Heatshield. 



s, in. 


(b) Aftbody and sting. 


Figure 6.40: Two-dimensional capsule surface heating after point insertion with cur 
vature clustering, solid=start, dashed=20% inserted, dash-dot-dot=benchmark. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.41: Two-dimensional capsule surface heating after point insertion with fluc- 
tuation minimization. solid=start,, dashed=20% inserted, dash-dot-dot=benchmark. 
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step of the adaption cycle. 

The curvature clustering cycle proceeded by deleting 1637 nodes, swapping 444 
edges, moving 696 nodes a total distance of 6.23 in., and adding 1599 nodes. The 
resulting heat transfer plots are shown in Figure 6.42. On the heatshield the stag- 
nation region heating has been improved and is now in excellent agreement with the 
benchmark. There is no further change in heating over much of the forebody between 
s = 0.2-0. 9, but the heating spike at the shoulder is more tightly resolved. In the 
wake region there is little change in heating on the aftbody, while the sting heating 
is increased, trending aw T ay from the benchmark. 

The fluctuation minimization adaption successfully removed 1601 nodes, but ran 
into trouble again while swapping. Windside vortices were spawned in the stagnation 
region by an oscillating bow shock. In an effort to damp the solution the CFL 
number was reduced by an order of magnitude, without producing an improvement 
in the solution or eliminating the vortices. The number of edges to be swapped was 
reduced to | percent, but the same vortices and oscillating bow shock appeared. 
For this case edge swapping was omitted and the adaption cycle continued with the 
nodal displacement step, where 706 nodes were moved a total distance of 0.7 in. 
Finally, 1638 nodes were added, yielding the results of Figure 6.43. The forebody 
heating is generally improved toward the benchmark solution, although there is a 
high-frequency oscillation in the data starting at s = 0.6. Heating in the wake is only 
slightly changed, though the aftbody heating is improved to match the benchmark 
between s = 1.2-1. 8. 


Axisymmetric 

The axisymmetric adaption component testing starts from the the baseline solution 
and mesh, which contains 11,125 nodes. The intention is to remove nodes without 
altering the surface heat transfer rates. About 10 percent of the nodes were success- 
fully deleted without impacting the solutions using both curvature clustering, 1029 
nodes, and fluctuation minimization, 1079 nodes. A close-up of the fluctuation split- 
ting mesh is shown in Figure 6.44, where it can be seen that most of the nodes are 
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(a) Heatshield. 
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(b) Aftbody and sting. 


Figure 6.42: Two-dimensional capsule surface heating after full adaption cycle with 
curvature clustering, solid=start, dashed=adapted, dash-dot-dot=benchmark. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.43: Two-dimensional capsule surface heating after full adaption cycle with 
fluctuation minimization, solid=start, dashed=adapted, dash-dot-dot=benchmark. 
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Figure 6.44: Coarsened axisymmetric mesh using fluctuation minimization. 


removed from the freestream while a smaller amount are removed from the wake re- 
gion. Further point reductions were tried, and the curvature clustering successfully 
removed another 1136 nodes with minimal change in the surface heating. However, 
removing 1081 more nodes with fluctuation minimization leads to the problem of a 
loss of shock capture, as was seen in the two-dimensional results of Figure 6.33. In 
this case, the solution deteriorates to an unacceptable level because of the bow shock 
blowout. 

Edge swapping is started from a mesh four times coarser than the baseline, tri- 
angulated from 63 x 45 = 2835 nodes, looking for the heating predictions to improve 
toward the baseline results. The mesh has 8290 edges, and 78, about 1 percent, 
with the largest error estimates were swapped using curvature clustering. Figure 6.45 
shows the initial and adapted heating rates compared with the baseline solution. The 
edge swapping produces only marginal changes in the heating rates. Fluctuation 
minimization swapped 88 edges and also show's minor changes in the heating, Fig- 
ure 6.46, except at the stagnation point where the adapted heating is dramatically 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.45: Axisymmetric capsule surface heating after edge swapping with curvature 
clustering, solid=start, dashed=adapt,ed, dash-dot-dot=DMFDSFV baseline. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.46: Axisymmetric capsule surface heating after edge swapping with fluctua- 
tion minimization, solid=start, dashed=adapted, dash-dot-dot=flnctuation splitting 
baseline. 
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Figure 6.47: Shock kinked in toward stagnation point at the :r-axis after axisymmetric 
edge swapping with fluctuation minimization, velocity contours. 

raised well beyond the correct result. Looking at the bow shock in the vicinity of the 
nose, Figure 6.47, a sharp kink in toward the body is seen along the axisymmetric 
axis. This kink in the bow shock elevates the pressure along the axis and causes the 
significant stagnation point heating over-prediction. 

The nodal displacement tests start from the same 2835-node meshes and solutions 
that were used to start the edge swapping check. Five percent of the nodes were 
moved using each scheme without producing a noticeable change in the solutions. 
A second pass of node movement was made, this time moving 10 percent of the 
nodes. Curvature clustering moved the nodes a total distance of 2.88 in. with an 
RMS of 0.017 in. The movement was performed at the bow shock and in the forebody 
boundary layer. The results, shown in Figure 6.48, show minimal change with a slight 
smoothing of the prediction on the heatshield between s = 0.2-0. 5. With fluctuation 
minimization the nodes were moved a total distance of 0.18 in. with an RMS of 
0.0015 in. Most of this movement was in the shock near the axis. Figure 6.49 shows 
basically no change to the heating levels due to the nodal displacements. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.48: Axisymmetric capsule surface heating after nodal displacements with 
curvature clustering, solid=start, dashed=10% moved, dash-dot-dot=DMFDSFV 
baseline. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.49: Axisymmetric capsule surface heating after nodal displacements with 
fluctuation minimization, solid=start, dashed=10% moved, dash-dot-dot=fiuctuation 
splitting baseline. 
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Point insertion is applied to the coarse 2835-node mesh using both schemes to add 
10 percent more nodes. Curvature clustering added 283 nodes at the shock and in the 
windside boundary layer, but does not change the surface heating appreciably, Fig- 
ure 6.50. Fluctuation minimization added 287 nodes, also in the shock and windside 
boundary layer. While the average heating levels on the forebody remain roughly the 
same, Figure 6.51 shows a wild oscillation in the surface heating has been introduced 
with the additional nodes. The sting heating shows a small reduction, but does not 
exhibit a clear improvement. 

For the full axisymmetric adaption procedure a cycle is defined similar to the 
cycle for two dimensions. Starting from a converged solution on a mesh with half the 
number of nodes as the baseline, 125 x 45 versus 125 x 89, 10 percent of the nodes are 
deleted, 1 percent of the edges are swapped, 5 percent of the nodes are moved, and 
then 10 percent more nodes are added back. The solution is re-converged between 
each step of adaption. 

Curvature clustering deleted 560 of the 5625 nodes, primarily from the freestream 
but also from the wake and in the sting/vehicle juncture. Of the 14,941 remaining 
edges, the 150 with the largest error estimates were swapped. These edges w r ere 
located in the bow shock and in the forebody boundary layer. Then the 253 nodes 
with the largest error estimates were moved. These nodes were also in the shock 
and forebody boundary layer. The nodes were moved a total distance of 1.2 in. with 
an RMS of 0.007 in. Finally, 570 nodes were added, enriching both the shock and 
the boundary layer. The results of the adaption cycle are shown in Figure 6.52. 
Essentially no change in surface heating is seen due to the adaption, other than a 
smoothing of the originally oscillatory data on the sting. 

Fluctuation minimization ran into problems with deleting nodes from the starting 
mesh for the full adaption test. Deleting 10 percent of the nodes proved to be too 
aggressive, leading to a loss of the bow shock capture. So the adaption cycle was mod- 
ified to only delete and insert 5 percent of the nodes, rather than the original target of 
10 percent. Deletion was successful removing 263 points entirely from the freestream. 
Edge swapping was viable for this case as 159 edges were swapped. The problems en- 
countered during the two-dimensional edge swapping were not encountered with the 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.50: Axisymmetric capsule surface heating after point insertion with curva- 
ture clustering, solid=start, dashed=10% inserted, dash-dot-dot=DMFDSFV base- 
line. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.51: Axisymmetric capsule surface heating after point insertion with fluc- 
tuation minimization, solid=start, dashed=10% inserted, dash-dot-dot=fluctuation 
splitting baseline. 
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(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.52: Axisymmetric capsule surface heating after full adaption cycle with cur- 
vature clustering, solid=start, dashed=adapted, dash-dot-dot=DMFDSFV baseline. 
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axisymmetric case. The nodal displacements were performed primarily at the shock 
with a small amount moved in the boundary layer, for a total of 268 nodes moved a 
distance of 0.22 in., RMS of 0.002 in. Point insertion added 280 nodes, mainly at the 
shock but also a small amount in the forebody boundary layer. The results of the 
adaption cycle are shown in Figure 6.53. Most of the heatshield remains unchanged 
except for an improvement between s = 0.3-0. 7 and a change in the stagnation point 
trend from being a large under-prediction to being a large over-prediction. Heating 
in the wake is largely unchanged other than an increase at the tail end of the sting. 


6.5 Discussion 

The primary new development of this chapter is the extension of the characteristic 
alignment mesh adaption strategy to the Navier-Stokes system of equations by a 
generalization of the fluctuation minimization formulas originally developed for the 
scalar advection-diffusion equation. Particular emphasis is placed on formulating the 
nodal displacement forcing functions for the upwind, positive, linearity preserving 
fluctuation splitting distribution scheme consistently in both axisymmetric and two- 
dimensional coordinates. 

The main test is a hypersonic, perfect gas wind tunnel case for a Martian probe 
model, configured for both the real axisymmetic geometry and the corresponding 
two-dimensional geometry. In considering the unadapted results, the two-dimension- 
al baseline solutions for both fluctuation splitting and DMFDSFV are obtained on the 
125 x 257 mesh, the same mesh for which the LAURA benchmark is grid converged, 
aside from at the stagnation point. The fluctuation splitting surface pressures are 
in excellent agreement with the benchmark while the DMFDSFV surface pressures 
match the benchmark on the heatshield but do not match well in the wake. For the 
forebodv heating fluctuation splitting shows excellent agreement with the benchmark 
aside from at the stagnation point while the DMFDSFV result is consistently high. 
The aft heating is about the same for both schemes, reasonably agreeing with the 
benchmark. 



6.5. DISCUSSION 


221 



(a) Heatshield. 



(b) Aftbody and sting. 


Figure 6.53: Axisymmetric capsule surface heating after full adaption cycle with fluc- 
tuation minimization, solid=start, dashed =adapted, dash-dot-dot=fiuctuation split- 
ting baseline. 



222 


CHAPTER 6. SYSTEM MESH ADAPTION 


For the axisymmetric unadapted results, the baseline solutions for both unstruc- 
tured schemes are obtained on the 125 x 89 mesh, which is a factor of four coarser 
than the mesh that Hollis had to use to get good results with the structured finite 
volume code NEQ2D. DMFDSFV surface pressures are excellent on the forebody, 
very good on the aftbody, but a little high on the sting between s = 2. 5-3. 5. The 
fluctuation splitting pressures are a little low on the f’orebody and excellent in the 
wake, splitting the difference between the two benchmark datasets. The fluctuation 
splitting forebody heating is an excellent match with the benchmark, aside from at 
the stagnation point. DMFDSFV over-predicts the forebody heating. Both schemes 
produce an excellent match with the benchmark in the wake, splitting the difference 
between the two computational benchmarks on the sting. 

In testing the component adaption techniques point deletion is the most beneficial 
tool. Coarsening works well removing up to about 20 percent of the nodes, but when 
taken too far leads to a loss of the bow shock capture, and a rapidly deteriorating 
solution. The other adaption techniques generally have minimal or negative impacts 
on the solution, with the exception of point insertion with fluctuation minimization 
where a little solution improvement occurs. Curvature clustering produces nodal dis- 
placements on the order of ten times larger than fluctuation minimization for the 
same number of points moved, but the fluctuation splitting solutions on the con- 
torted meshes are less robust. The lack of robustness is usually caused by distortions 
to the bow shock near the symmetry axis. Common failure modes are an oscillating 
shock shedding windside vortices or a steady shock kink producing unexpected pres- 
sure contours and surface heating spikes. The root cause is that the local, compact 
adaption techniques employed here are mis-aligning the mesh with the bow shock. 
The typical behavior is that the fluctuation splitting scheme is more sensitive to grid 
perturbations introduced by the fluctuation minimization adaption in corrupting the 
solution, but is generally more accurate than DMFDSFV on suitable meshes. The 
component tests indicate that the mesh adaption for systems can not be nearly as 
aggressive as was possible for scalar equations. 

For the full adaption cycle in two dimensions with curvature clustering the heating 
in the stagnation region and at the shoulder improves, the sting heating is worse, and 
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the rest of the vehicle sees no change. With fluctuation minimization the forebody 
heating improves but becomes oscillatory, the aftbody is the same, and the sting 
improves. For the axisymmetric case the curvature clustering cycle does not change 
the magnitude of the heating but does smooth the initially rough heating on the sting. 
Fluctuation minimization produces some improvement over a quarter of the heatshield 
but also makes the sting heating a little worse at the tail end, leaving the rest of 
the heating unchanged. Although heat transfer is a variable of primary interest for 
aerothermodynamics and may be a first consideration for driving the adaption criteria, 
the present results emphasize that for blunt-body hypersonic aerothermodynamics 
the bow shock placement is critical to the forebody heat transfer rates. Due to the 
strength of the hypersonic bow shock, discrete jumps in shock position from one node 
to another produce significant pressure disturbance waves which strongly affect the 
flowfield in the subsonic bubble and can overwhelm more subtle adaption within the 
boundary layer. 

To conclude, the adaption techniques do not significantly improve the solutions, 
the adaption effectiveness is not consistent, and the solutions display a lack of robust- 
ness. Point deletion works the best, but is risky if applied too aggressively. These 
results indicate that the extension of characteristic alignment based adaption from 
scalar problems to the Navier-Stokes system is not straightforward due to the non- 
uniqueness of the characteristics and the complexity of the flowfield features which 
are not adequately modeled by the minimization of a single functional based upon 
fluctuations as formulated here. 

These local adaption strategies are designed to be synergistic with the compact 
solvers, and are attractive candidates for parallelization, but the present results may 
be indicating that the adaption requires a global view to introduce a level of smooth- 
ness. Dr. Peter Gnoffo[123] has suggested reworking the adaption strategies to in- 
crease the support of the discretization stencil. Such an increased stencil may allow 
for the introduction of smoothing parameters or other grid-quality controls, but at the 
expense of requiring more than nearest-neighbor communications, complicating any 
extension to a parallel computing strategy. The upper limit on an increased stencil 
would be a fully implicit treatment. 
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As a post-note to this chapter. Yamaleev and Carpenter[124] have recently pre- 
sented a very detailed mathematical analysis of the effect of grid adaption for super- 
sonic blunt bodies in the context of second and fourth-order structured finite volume 
schemes. They specifically consider nodal displacements and point insertion, conclud- 
ing, “The grid refinement study shows that for the second-order scheme, neither grid 
adaptation strategy improves the numerical solution accuracy. . . ” For the fourth- 
order scheme the situation is worse because, . . the design-order error component 
drastically increases because of the grid nonuniformity.” A consensus on adaption 
guidelines for general supersonic flows with captured multi-dimensional discontinu- 
ities remains elusive, for both structured and unstructured meshes. 



Chapter 7 

Summary and Recommendations 


This report addresses the suitability of fluctuation splitting for hypersonic aerother- 
modynamics. A desirable aerothermodynamic solver must be robust in the presence 
of strong shocks, meaning that the solution must remain stable and positive with- 
out radiating dispersion or phase errors. The desirable scheme must also be of high 
accuracy to resolve the derivative quantities heat transfer and shear stress, and thus 
should be non-diffusive. The classic approach to reconciling these two criteria is the 
construction of a non-linear second-order scheme with reduction to first-order in the 
vicinity of captured discontinuities. Perhaps the most popular of these schemes for 
unstructured meshes is finite volume with Roe upwind flux difference splitting and 
limited linear reconstruction (referred to as DMFDSFV here). A drawback to this 
scheme is the reliance upon a locally one-dimensional approximate Riemann solver 
even when applied in two and three dimensions. Fluctuation splitting possesses the 
benefits of the traditional scheme and adds a truly multi-dimensional upwind distri- 
bution. But the discussion of a flow solver independent of the computational domain 
is only of academic interest, so this report also seeks to develop a mesh adaption 
scheme in conjunction with the desirable solver to improve the results and reduce the 
expense of analyzing a vehicle geometry. 

The fluctuation splitting scheme considered here had been introduced by Sidilkover 
as the optimal compact zero-cross-diffusion solver for linear hyperbolic equations, 
with an extension to the two-dimensional Euler equations on Cartesian grids. The 
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formulation lacked eigenvalue limiting, an axisymmetric derivation, and a viscous 
coupling, had not been demonstrated on a general unstructured mesh, was untested 
for non-linear advection-diffusion, had not been applied to a hypersonic problem, had 
not been used to calculate heat transfer, and had never been evaluated relative to the 
current most popular strategy — DMFDSFV. The evaluation of fluctuation splitting 
begins by proving the equivalence of DMFDSFV and fluctuation splitting in one 
dimension along with a strategy for viscous coupling to the inviscid flux distribution. 
Then the extension of both schemes to scalar two-dimensional problems is detailed to 
underscore the fundamental differences between a multi-dimensional upwinding and 
the locally one-dimensional treatment, with particular emphasis on the production of 
artificial dissipation. The superiority of fluctuation splitting over DMFDSFV" for both 
linear and non-linear advection-diffusion is demonstrated and a new mesh adaption 
strategy for scalar problems is developed to exploit characteristic alignment with 
fluctuation splitting. This adaption scheme maintains the fully local, compact stencil 
of fluctuation splitting to allow for future parallelization and use on computer clusters, 
and is able to produce accuracy and efficiency gains for the scalar problems. 

The details of extending multi-dimensional upwind fluctuation splitting to a prac- 
tical aerothermodynamic solver for the two-dimensional and axisymmetric Navier- 
Stokes equations are detailed and a suite of verification and validation cases compares 
the fluctuation splitting and DMFDSFV" schemes. Fluctuation splitting is still to be 
preferred over DMFDSFV for the system of equations, more so for viscous problems 
than inviscid, but the difference is not as dramatic as is seen for the scalar model 
problems. While fluctuation splitting should certainly be considered as a viable alter- 
native to DMFDSFV when selecting algorithms for code development, the differences 
in capability demonstrated here are not significant enough to obsolete existing codes 
based upon DMFDSFV" technology. 

The characteristic alignment mesh adaption strategy is extended to two-dimen- 
sional and axisymmetric systems as a generalized form of fluctuation minimization. 
The practical implementation for systems proves to be very difficult, and in the cur- 
rent adaption framework only the point deletion strategy is recommended. None of 
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the other adaption components produce significant accuracy improvements and of- 
ten lead to a loss of robustness of the solver. For inviscid or subsonic flows where 
heating is not an issue, perhaps the current adaption could be of use. For aerother- 
modynamics though, the local adaption strategy as it is currently presented is not 
recommended. 

Parallelization of the fluctuation splitting scheme appears to be attractive because 
of the compact stencil. The weak implementation of the boundary conditions devel- 
oped here should be an enabling mechanism for exchanging boundary data between 
partitioned sub-domains in a parallelized scheme. Another significant operational im- 
provement could be made by casting the fluctuation splitting in an implicit scheme, 
with either a point-implicit or colored strategy the most likely format. Two-equation 
field models for turbulence appear to be the best path for extending the current fluc- 
tuation splitting formulation to Reynolds and Farve-averaged flows. One remaining 
significant unknown for fluctuation splitting with regards to aerothermodynamic ca- 
pability is the inclusion of thermal and chemical non-equilibrium. In particular, the 
transformation of the flux Jacobian to auxiliary variables may not yield an easily 
f'actorizable system. Appendix B expresses the linearizations and transformations in 
terms of a general gas, but all the details of such an extension need to be worked out 
before a fair assessment of implementing a non-equilibrium scheme can be made. 

Modifications to the adaption strategy for increased robustness should relax the 
constraint that the adaption be as compact as the solver, focusing on some form of 
increased stencil as a means to introduce smoothing. Perhaps, though, a completely 
different adaption strategy is called for, something with a scaffolding analogy by which 
the bow shock, embedded shocks, vortices, and other inviscid phenomenon are framed 
with lines or shell elements while the gaps in the mesh framework are filled in with 
an unstructured triangulation. Whatever adaption strategy is pursued, the present 
results underscore the need for excellent bow shock alignment as a prerequisite to 
boundary layer refinement for the hypersonic blunt-body problem. 
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Appendix A 
Limiters 


A limiter is a function used to limit the ratio of two values, satisfying, 


V>(0) =0, ?/’(!) = 1 

A symmetric limiter is defined by, 



(A.l) 


(A.2) 


Symmetric limiters can also be expressed in terms of symmetric averaging functions, 
M^, obeying, 

q ip (~) = M it {p,q) = Mip(q,p) =pip(^J (A.3) 

A limiter that can achieve a value greater than unity is termed a compressive limiter. 


Degenerate limiters 

Two degenerate, non-symmetric limiters, satisfying only one of the constraints in 
Eqn. A.l, are useful. The first order limiter, ib = = 0, is employed to limit a 

scheme to first order spatial accuracy. In contrast, an unlimited scheme results from 
the choice ?/; = 1 . 
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Minmod 

The Minmod limiter[125] is a non-compressive, symmetric limiter defined as, 

f N 

P 


or, 


ip ( - ) = ma.x(0, min(l, p/q )) 


(A.4) 


The associated averaging function is, 
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The Minmod limiter is the non-compressive limit of a generalized limiter of 
Sweby[126]. Minmod is achieved by setting the parameter e = 1. The upper limit on 

e is the Superbee limiter[127], e = 2. 

/ \ 

P \ rn i \ _\1 (A 7) 


')/’(-) = max[0, min (ep/q, 1), min (p/q, e)] 


, , P 

ip 


0 

pq< 0 

zp/q 

e\P 1 <\q\ 

1 i 

f \p\ <kl< e\ p 

p/q 

kl< \p\ < e\ q 

£ 

e|?l< bl 

0 

pq < 0 

zp 

£\p\<\q\ 

q if 

\p\<\q\< £\p\ 

P 

\q\ < \p 1 < e|</l 

zq 

£\q\ < \p\ 


(A.8) 


M (p, q) = { 


Similar, non-symmetric limiters have been proposed by Chakravarthy[128], 

f P 


(A.9) 


ip 


= max[0, min (p/q, e)] 


(A.10) 
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and Barth [22], 


4> 


= max[0, min (sp/q, 1)] 


(A.ll) 


The upper bound on these limiters is 1 < e < 2. 


Van Leer 

The harmonic Van Leer limiter [4] is a symmetric compressive limiter with an upper 
bound of 2. 


, , P 


f+lfl 

i+ISI 




= pq+\pq I . 

q 2 +\pq\ 

\p\q + p\q I 

bl+M 


bl+klf 

bl+M 


(A. 12) 
(A. 13) 


Van Albada 


The Van Albada limiter[129] is a symmetric, differentiable compressive limiter with 
an upper bound of 1.18. 


</>i - 


p 2 +g 2 i v 

q'-—c~ q 
p 2 +s 2 


1 + 


M, 


q 2 +e l 


p 2 + £ 2 + ( q 2 + £ 2 )\ 
p 2 + £ 2 + q 2 + £ 2 


(pq + £ 2 ){p + q) 


(A.14) 


(A.15) 


p 2 + q 2 + 2£ 2 

For this limiter the small parameter, e, serves to reduce the limiting in smooth 
regions. Assuming 0(1) variations over a normalized distance, this parameter is 
scaled as £ 2 ~ a 3 , where ^ is the mesh spacing. 

Van Albada presents the averaging function as in Eqn. A.15. Limiter-function 
forms have been presented corresponding to this averaging function in popular sources, 
such Hirsch[125]. Denoting the ratios p/q by r and e/q by s, the limiter presented in 
these sources is, 


4>(r) = 


(r + l)r 
r 2 + 1 


which, lacking e, differs from Eqn. A.14, 


Hr) = 


(r + l)(r + -s 2 ) 
r 2 + 1 + 2.s 2 
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The more popular forms of this limiter do not turn off in smooth regions, thus missing 
an essential feature of the original Van Albada formulation. 



Appendix B 

Governing Equations 


B.l Compressible Continuum Gas Dynamics 

The Navier[87]-Stokes[88] equations of mass, momentum, and energy conservation for 
compressible gas dynamics are presented. The formulations are in terms of Cartesian 
coordinates, but include an axisymmetric source term, 1 * B, which is set to zero for 
two-dimensional and three-dimensional applications. 

The derivation of the equations assume a Newtonian stress-strain relationship, 
Stokes hypothesis on the bulk viscosity, Fourier’s law for heat conduction, no body 
forces, and no external heat addition. 

Non-dimensionalization is performed for the following quantities as: length, L re j. 
velocity, T*,, time, L ref /V t x , energy, V£, density, p^, pressure, pooV^, viscosity, p^, 
temperature, T, thermal conductivity, p^V^/T^, and gas constant and specific 
heats, V^/Too. 

The Reynolds number appearing in the non-dimensional equation is defined, 

R eoo = p °° V °° L ™f (B.l) 

Poo 

1 The axisymmetric source terms were derived from the orthogonal curvilinear formulations of 

Back[130] and Anderson, Tannehill, and Pletcher[89]. 
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B.1.1 Vector Notation 

The system can be expressed in conserved variables as, 

Uj + V-F =-wB (B.2) 

where F = F { — F u , B = B* — B u , and, 
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B' = B' - B' 
B'‘ = (0, 0, P, 0) T 

B'" = -=-?- (o,0,— 

R,.„ V y 3 


(B.ll) 

(B.12) 

(B. 13) 


Also, when using the form Eqn. B.9 the divergence of the velocity in the viscous 
terms, and B' , is replaced by, 


- T -> 1 ^ VN7zu a - T -t vjv 

V-V — » — V-(w a V) = V - \ H = V V H 

^ a t&a V 


(B.14) 


B.l. 2 Two-Dimensional/ Axisymmetric Indicial Notation 


The axisymmetric source terms are included in braces {}. 

Continuity: 

Pt + (PU) X + (pV)y =~ {/*>} 

TXaPt + ZUa{pu)x + &a{pv)y = 

w a p t + zu a {pu) x + {w a pv)y =0 

Momentum-x: 

(. PU)t + ( PU 2 + P) X + ( PUV)y = 0[/i(2U X - Vy)] x + [fi(Uy + V X ) ] y ^ 

w ( i r / . 2 

+ — i -puv + p\u y + v x ) - -{pv) 3 


(B.15) 

(B.16) 

(B.17) 


(B.18) 


vo (2 

&a(pu) t + W a {pU 2 + P) x + VJ a{puv) y = ~[//(2 U x - V y )\ x + [p{u y + W x )]j 


+ G7 | -puv + -2— p(u y + v x ) - \{pv) x | (B.19) 


TXa{pu)t + ™ a {pU 2 + P)x + ( WaPUv\ 


Jj— f ^a[p{2u x - Vy)\ x + [m a p(u y + Ux)]^ - ^-{u v h (B.20) 
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Momentum-y: 


(pv)t + ( PUV) X + ( PV 2 + P)y = 


R P 


[p{Uy + ] a; + Vy U x )\y 


W [ 9 1 

+ - - P V + — 

V { Re a 


O 2 ( \ 4 P V 

2 P v y - 3 (Hv - 3 ~ 


(B.21) 


&a(pv)t + W a (puv) x + OT a (pW 2 + P) y = ([p(u y + V x )] x + ^[p(2v y - U x )\ 


+ VJ < —pv + 


1 


R e 


2 4 pv 

2pVy - -{pv)y — 


(B.22) 


'UD ( 2 

ZU a (pv)t + ZU a (pUv) x + (W a pV 2 )y + W a P y = [p{u y + V x )\ x + ^[p(2v y - U x )] j 

Re<x V 4 


+ 


2ro 


3 R, 


’Coo 


Energy: 


(p£)t + ( puH) x + (pvH) y = 


R , 


Coo 


{kT X ) X + {nTy)y + 


— 2— j (B.23) 


P ( ~u{2u x - Vy) + V^y + V x ) 


+ 


w \ 1 

+ — \ —pvH + — 

y { R 


P j U{Uy + V X ) + ~v(2 Vy - U x ) 


(B.24) 


UTy + p(u(u y + V x ) + 2 VVy) 

'Coo 

2 4 

- ~{u(pv) X + v(pv)y) - ~pv(u x + Vy) 


VJ 

™ a{pE)t + W a (puH) x + W a (pvH)y = 

ilp 


{kT x ) x + (kTj,)j 


+ 


P ( 3^(2 u x - Vy) + v(u y + v x ) 


P ( U(Uy + V x ) + -v{2 Vy - U x ) 


+ 


Pm < -pvH- 


J y 


(B.25) 


R, 


2 

- -(u(pv) X + v(pv)y) - -pv(u x + Vy) 


nT y + p,(u(u y + v x ) + 2vv y ) 
4 
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a(pE) t + ru a (puH) x + (■ w a pvH) y = —— w a {nT x ) x + (w a K,T y ) y 

ilp 


Tw a p ( -u(2u x - Vy) + v(u y + v x ) 


+ zu a p I u(u y + v x ) + -v{2 Vy - u x ) 


{{puv) x + (pv 2 )y} 


(B.26) 


B.2 Linearizations 


The inviscid flux Jacobian is defined as, 


(B.27) 


For a general gas in two dimensions with, 

op _ op 

dpu U dpE ' 

the Jacobian is, 


dP_ _ _ dP 
dpv 1 dpE 


(B.28) 


A x = 


P p - v 2 u( 2 - P pE ) -vP p 


(B.29) 


u(Pp - H) H - V 2 Pp E -uvPpE u( 1 + P pE )_ 


A y = 


Pp - V 2 -uPpE v{2 - PpE) PpE 

_v(P p — H) UV Pp E H - V 2 PpE v{l + Pp E )_ 


(B.30) 


The projection of the Jacobian onto a unit vector can be written in compact notation, 


A -h = P p h T -VV T VI + V T h - n' VPpE h T P pE 
V(P, -H) Hn VVP pE V(1 + ) 


(B.31) 
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where the projected velocity is V = V -h. 

Specializing to an ideal gas with an equation of state, 

{puf + (pv) 2 


P = (7 - 1 )pe = (7 - 1 )[pE- 


2 P 


P '= 2 


1 1 y2 

v ) 


PpE = 7-1 


where V 2 = V -V = u 2 + v 2 , the Jacobian becomes, 


A x = 


- u 2 


1 0 0 

—u( 7 — 3) —v(j — 1) 7 — 1 

v u 0 

uQy-V 2 — H) H — u 2 ( 7 — 1) — uu( 7 — 1) 7U 


-uv 


A y 


—uv 

^V ' 2 - u 2 


0 

v 

-u(7 - 1) 


1 

« 


0 

0 


-v(7 — 3) 7 — 1 


u ( I 2^^ r2 — H) —uv(j — 1) H — u 2 (7 — 1) 7« 


The projected Jacobian is, 


A -h 


0 n x n y 0 

— uV + 2 ^-V r2 n x V — (7 — 2) un* un y — (7 — l)tm x (7 — l)n x 

— t'V + 2 ^2-V' _2 n y un x — (7 — l)im y V — (7 — 2)un y (7 — l)n y 

Lv( 2 f i T' 2 - J H‘) ffn x - (7 - 1)«V Hn y — (7 — l)vV jV 


The eigensystem for the projected Jacobian is formed as, 

A-h = XAX -1 


The eigenvalues are, 


A = diag(V, V, V + a, V - a) 


(B.32) 

(B.33) 

(B.34) 

(B.35) 

(B.36) 

(B.37) 

(B.38) 


with the sound speed defined, 
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The associated eigenvectors are the columns of X, 

10 1 1 

u — n y u + an* u — an* 

(B.40) 

v n x v + an y v — an? 

■rr rm x — un? H + aV H — aV 

B.3 Variable Transformations 

Transformations of the dependent variables can be performed from the conserved 
variables, Eqn. B.3, to the primitive variables, 




satisfying dV = U v dV, dV = V, dXJ. and Uy 1 = V f/ . For a general gas, the 
transformation matrices are, 

1 0 0 " 

U v - = F T pi 0 

l ' 2 - pV -L- 

*f)E *pE _ 

1 0 0 

Vf 7 = -±V T i/ 0 

p p 

Pp - V f’pB 

While a perfect gas assumption leads to, 

1 0 0 

Uy = y T pi 0 (B.44) 

— pV — 

2 1 7-1. 

1 0 0 
V f/ = -iF T 0 

_^V 2 -( 7 -l)V 7 -l_ 


(B.42) 

(B.43) 


(B.45) 
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The projected flux Jacobian, Eqn. B.31, can be transformed as, 

h-Ap = V jj ^A-n') V v 

For a general gas, Eqn. B.46 leads to, 

V ph 0 

h-Ap =0 VI - p n T 

_0 p[P p - P pE (V 2 - H)]h V _ 

while the perfect gas version is. 

V ph 0 

h-Ap = 0 VI -h T 

p 

0 7 Ph V 

Introducing the auxiliary variables such that, 


(B.46) 


(B.47) 


(B.48) 



(B.49) 


where ds = dp—\dP , a transformation from the primitive variables can be introduced 
satisfying dV = Vyr-dW, dW = W v'dV, and Vj^ 1 = Wy, with, 

1 0 4 

a z 

Vw = 0 i/ 0 (B.50) 

0 0 1 

i o -4 

a z 

Wy = 0 pi 0 (B.51) 

0 0 1 
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The projected flux Jacobian can be transformed accordingly, 

h-A = W v {n-Ap)V w 

V {1 - ^[ Pp - PpE (V 2 - H)}} h 0 

= 0 VI n 1 

o [P p - P pE (V 2 - H)]n V 

V 0 0 

0 VI h T 
0 a 2 h V 

The eigensystem of Eqn. B.54 can be expressed as, 

n-A = XAX~ l 

where the eigenvalues remain the same as for Eqn. B.38, 

A = diag(V, V. V + a, V - a) 

and the eigenvectors are, 

10 0 0 

0 n y n x n x 
0 — h x n y n y 

0 0a —a 

1 0 0 0 ' 

0 h y —h x 0 

0 ^ 

0 I”" ¥? 2a- 

The conversion from conserved to auxiliary variables is, 

dU = Uy dV = Uy V w dW = \J W dW 
dW = Wy dV = W y V if dU = W„ dU 




(B.52) 

(B.53) 

(B.54) 

(B.55) 

(B.56) 

(B.57) 

(B.58) 

(B.59) 

(B.60) 
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For a perfect gas, these matrices are, 


W v 



1 

0 


1 

a 2 

= 

1- 

I 


a- 


V 2 

V 


1 To 


2 


7 — 1 T 

_ To 


V 


J 

r 


h 


-V T 


I 



2f±V ' 2 —(7 — 1)V 7-1 


where the enthalpy is obtained from the total enthalpy, 



and the stagnation temperature is defined, 


n 

T 


= 1 + 


7 ~ 1 
2 


M 2 


(B.61) 


(B.62) 


(B.63) 


(B.64) 



Appendix C 

Axisymmetric Source Term 
Integration 


C.l Inviscid 


The perfect-gas inviscid axisymmetric source term can be written either in the form 
of Eqn. B.7 as, 



(C.l) 


or in the alternate form of Eqn. B.12 as, 


B' = [0, 0, P, 0] T 


0 , 0 , (yZ\Z<± 




(C.2) 


The parameter vector, Z, is assumed to vary linearly over a triangular element il. 
The integration of B over i} is performed by subdividing into thirds along the 
median-dual mesh, as depicted in Figure C.l. Notice that each sub-element 12, is a 
quadrilateral over which the parameter vector varies linearly. Also, the area of each 
is one-third of the area of the triangular element, 



(C.3) 
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Figure C.l: Subdivision of triangular element into three quadrilateral integration 
areas. Dashed lines are the median-dual mesh. 


A linear variation, say of the m fh element of Z, can be represented as, 


Z m {x,y) 


2 (®m A b m x + c m y ) 

-^UjkZmj [(£ - Xi)(y k - yi) + (y - yi){xi - x k )] 


(C.4) 


where, 


dm = Z mi (x 2 y 3 - x 3 y 2 ) + Z m2 (x 3 y x - x 1 y 3 ) + Z m , i (x x y 2 - x 2 y\) 

— UjkZmj (x k yi Xiy k ) 


(C.5) 


b m = Z m i (y 2 - y 3 ) + Z m2 (y 3 - rq) + Z rni {y-y - y 2 ) 

£ijk.Z m j (| Ik £/?') 


(C.6) 


dm Z l!:i (.rq 3^2) 4 “ Zjji2 {,X\ X3) + Z ma (x 2 x\) 

dijkZmj {X'i iC/j) 


(C.7) 


In this appendix, eij k is used to represent the cyclic permutation summation operator. 
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C.1.1 Integration of Eqn. C.2 

The only non-zero term in Eqn. C.2 is B' 3 = P. The integration over a quadrilateral 
sub-element can be broken out as. 


■In, 


Pdi\ =- — - / Z\Z\ — 


/ Z X Z 4 di\ - ^ / Z\ + Z\ dili 

Jn, 


iQi 7 ./a, 27 J a . 

For a linear variation of Z, a typical term of Eqn. C.8 can be written, 


Jn, 


dili 


1 


J + b m' x + c m y)(a n + b n x + c n y) di\ 

I T (<l„il>n T (7 n b m )x T 0, n Cm ) V 

Jn. 


454 JQ Z 


(C.8) 


(C.9) 


T {b m c n T b n c m )xy T b m b n x T c m c n y dPli 


Equation C.3 leads directly to. 


L: 


/o. - 

(X777 ( 3 't? (Iu Lj Chr~n CL' 


(C.10) 


The remaining integrals in Eqn. C.9 are evaluated using a bilinear mapping to a 
unit square in (£, p) space. Introducing the notation i— and i+ to represent the node 
clockwise/counter-clockwise, respectively, from node i on triangle (l, ie, 


i— i i+ 

3 1 2 

1 2 3 

2 3 1 

the mapping from ( x , y) to (^, 77) is, 

x(^ rj) = 77 + r 2 ^ + r 3 rj + r^rj 

y(€, v) =h + t 2 £ + hji + t^rj 


where the coefficients are, 



1 

2 

3 

4 


r 

Xi 

+ 

1 

\ (Xi- ~ Xi) 

\{~Xi - + 2 Xi 

- x i+) 

t 

Vi 

\{Vi+ - Vi) 

\{Vi- ~ Vi) 

\{-Vi- + 2 yi 

~ Vi+) 


(C.ll) 

(C.12) 
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£ runs from node i, £ = 0, to the midpoint of side ii+, £ = 1. ?/ runs from node i, 
r] = 0, to the midpoint of side ii—, ij = 1. (£, rj) = (1, 1) corresponds to the centroid 
of il. The inverse Jacobian of the transformation is, 


J 1 =x^y Tt -x v y^ 


St 

2 


1 - 1 (f + v) 


and a transformed differential area element is, 


(C.13) 


di\ = J 1 dr) 


(C.14) 


Using the triangle centroid coordinates, 

Xi + X 2 + :l'3 


X = 


y = 


y\ + V2 + 2/3 


3 " 3 

the integrations in Eqn. C.9 take the following forms, 


xdili = + 5xi) 

36 


JVi 


X ydi\ = ^ 


Jn, 

f Sj 

/ ydtoi =^(7y + 5yi) 

Jn, 36 

3 

Uxy + 11 (x t y + xyi) + 9 xffi + ^ Xjyj 


3 = 1 


x 2 dfij = — 1 - ( 14x 2 + 22a:jf + 9 .t 2 + V" 

1 zlzl V ' 


y 2 dili = 




144 

5V 

144 


Xi 


3= 1 
3 


14 y 2 + 22 y(y + 9 y 2 + ^ ; 


i=i 


(C.15) 

(C.16) 

(C.17) 

(C.18) 

(C.19) 

(C.20) 


C.1.2 Integration of Eqn. C.l 


The integration of the axisymmetric source term, as given in Eqn. C.l, over a quadri- 
lateral sub-element is, 



/ ~ T3 Z m dl\ 

Jn i V 

7^2 / _ (°3 + + oy) (a m + b m x + c m y ) dil, 

Jn. ; V 

BTi -I- BT-2 + BT ?J + BTi 


(C.21) 
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where, 


BTi 

/ 

45| A 

dili 

h 

(C.22) 

bt 2 

b 3 Cm 4~ b m c 3 I 

4>4j 

x di\ 

(C.23) 

bt 3 

c 3 c m / 

= y a ,» dn ‘ 


(C.24) 


IS? L, v 


+ (« 3 b rn + a m b 3 ) x 4- b 3 b m x 2 ] di\ 


Equation C.3 leads directly to the first integral, 

7 jrp f 4 “ 3 

tttt; 


(C.25) 


(C.26) 


Using the same bilinear mapping to (£,, )]) space as was described in the previous 
section, Eqn. C.23 can be written. 


BT 2 = ' j J in + r 2 $ + r 3V + nin) (j - ^ 

= &3C ^+^ C - 3 [24rj + 11 (r- 2 + r 3 ) + 5r 4 ] 

Similarly, Eqn. C.24 becomes, 

BT 3 = 9 gg^ T [24fi 4- 1 1 (h + t 3 ) + 5 f 4 ] 
BT 4 in Eqn. C.25 is further simplified, 

BT a = BT41 4 - BT42 + -BZ43 

where, 

nm d 3 Qjm f f //• 


b^Cm + b m c 3 f 


Z + V\ St 


~r d£ drj 


(C.27) 


(C.28) 


(C.29) 


BT 41 = 
bt 42 = 

BT 43 = 


I n? I u ' “l 

4Sf -hk y 

^ 4,,, 4- ci m b 3 I x 




b 3 b m /' x : 


/ — dJ] 
7m y 
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Using the coordinate mapping in Eqns. C.11-C.12 to transform Eqn. C.30 leads 


to, 


BT U = 


a^cim 

24S t 

245V 


ff 
I ln 

Jo 


3 — ^ — r? 


U + h( + hv + Uhl 

U + h + (£3 + £ 4 ) 7 / 


d£ drj 


U + hv 


1 


t 2 + t 4 i} 


3 - i] + 


1 


U + h'V 
t-2 + UV 


drj 


ln 

U 


1+t i 


(C.33) 


Options for evaluating the remaining integral in Eqn. C.33 are to expand the natural 
logarithm in a series expansion as, 


ln 


ti + h + (V3 + U)v 
t\ + hv 

OO 1 


= 111 


U + to 
t-l 


n =1 


2 n 


1 

(^3 + *4)77 

2n— 1 

hV 

2n— 1 'j 

l 

, 2 (ti + t 2 ) + (h + f 4 )d. 


. 2 + + t 3 r/_ 

) 


(C.34) 


or to resort to a numerical integration scheme, such as Gaussian quadrature. 
Equation C.31 can be transformed as, 


BT, 


42 


bm + (JmJ ) 3 

245V 

®3 bm + o, m b 3 
245V 


1 r'i + + r 3 r/ + r^y 

t\ + h€ + t-.Ji + 1 1 £// 

1 


./o ./O (*1 + hV) + (0 + Ud/K 


(3 — — 77 ) d£ drj 


jo Jo 
1 


(C.35) 


• {(n + 7 3 r/)(3 - y) + [( r 2 + 7 4 r/)(3 - ?/) - (ri + r 3 rj)] £ - (r 2 + 7 4 r/)£ 2 } d£ dr/ 
The components of Eqn. C.35 are evaluated as, 

V'i + 7*377) (3 - ?/) 


= 1,1 + , ' :!,,l . (3 ~ [,) In ( + 1 


Jo ( U + hv) + (h + Uv)£, ' h + Uv \ti + hV 

1 f 1 (7 2 + 7 4 r/)(3 - r/) - (rq + r 3 r/) 


(C.36) 


J 0 Jo 


(ti + f 3 7/) + (t 2 + 


d? dr] 


3r 2 - r ! 

u 


ln 

7 4^2 


f 4 




+ 1 


+ ( 37 4 - 7 2 - 7 3 ) 


I-4h 

t 4 


t 2 

C 4 


t 4 


+ 1 


-H + 

2f 4 t: 

r I 


2 In |t 2 + t 4 | — - 7 ^ 


(C.37) 


ln |t 2 + i 4 r/| dr/ 


[(7 2 + 7 4 r/) (3 - r/) - (n + 7 3 r/)] (U, + hv) ^ / 0 + Uv 


(t 2 + Uv ) 2 


U + hV 


+ 1 dr/ 
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r -2 + Ti T] 


e<K<hi 


o (t\ + hv) + (h + Uv)£ 

T 2 + r 4 ?? ( 1 _ tg + 6 3 7/ 

o * 2 + Uv l 2 t 2 + hn [ 


1 il + t 3 ?] ^1 + t -2 + (t 3 + ti)lj 


= — In 


2t 4 


l + 1 


r 4 r 4 t 2 , 

4 7T In 

2 U 2 t\ 


t 2 + t 4 r] 

+ 


t 1 


r 2 ti 


h + hi] 

r 2 h 


dij 


r 2 t 3 + r A ti 

t' 2 

l 4 


In 


u 

h 


+ 1 


+ 


t 2 


t 2 + 1 4 


-1 


t 4 (t 2 + £4) 

Tit 3 
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Integrating Eqns. C. 40 -C .43 with respect to £ and then ij, where possible, yields, 
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• W. A. Wood and W. L. Kleb, 2-D/ Axisymmetric Formulation of Multi-Dimen- 
sional Upwind Scheme. AIAA Paper 2001-2630, June 2001. 

Future publications are anticipated for the following topics: the placement of 
eigenvalue limiting within the fluctuation splitting context from chapters 3 and 5; 
the weak formulation of the boundary condition for the upwind fluctuation splitting 
distribution scheme from chapters 3 and 5; and the adapted and unadapted capsule 
results of chapter 6. 



Appendix E 

Navier-Stokes Solver and Mesh 
Adaption Code 


The source code, written in C, for the two-dimensional and axisymmetric fluctuation 
splitting and DMFDSFV Navier-Stokes solvers along with the curvature clustering 
and fluctuation minimization mesh adaption schemes is included. The codes are 
actually written as the text documents shown here, and are converted to executable 
code by a Perl script as part of a Makefile compilation. 
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raesh_adapt = none; Perform loops over the domain to evolve the solution in pseudo- time. 
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bottom_.ce 11 = ptr_new ; F0R(i,0,2){ 
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void kill_node( Node *node ){ 1.5.9 Number Nodes in a Stack 
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void hook_from_node_to_edge ( Node *node, Edge *edge 



if ( (*n2e = (Node2Edge *) mallocC sizeof (Node2Edge) )) ){ 
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void build_bedges ( Node *n0, Node *nl , int safe 



bedge->node [0] 
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F0R(i,0,4) printf ("V%d=% .5e, ", i, n->vars .primitive [i] ); putsC"" 
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Hand-built BLAS routines 
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Ascent of Man ++cin 
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node->links . node .number, sigm[l] ) ; length and the projected wavespeed along the edge. Also check the cell midpoint to handle 
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Characteristic Alignment Adaption 
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pia *= sqrt (cell_right->ge om. metric . area) ; 

cell_f luc (cell_right->geom . connect .nodes , fluctuation) ; 

edge_f luctuation += dnrm2(NEQS, fluctuation, 1) /pia; 



FOR C j ,0,NEQS) AminusB[i] [j] = Lower[i][j] = Upper[i][j] 
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dbetadx [3] [2] = a2 * (dbetadx[2] [3] = 1 /2); for( sum=0, j=0; j<i— 1 ; ++j ) sum += Lower[i][j] *y_is_Ux[j]; 
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u = celln->cell->tilde .primitive [1] ; 
v = celln->cell->tilde .primitive [2] ; 

v solve for ^ and ^ using forward and backward substitution. a2 = SQUARE(celln->cell->tilde.a) ; 



,_z[i] = node2->vars .parameter [i] - nodel->vars .parameter [i] ; Wzdzdy[i] = ddot (NEQS , WZ[i],l, dZdy,l); 
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dvdx [0] = (dZdx [1] -Z2[l] *dZdx[0] /Z2[0]) /Z2[0] 

dvdx [1] = (dZdx [2] -Z2[2] *dZdx[0] /Z2[0]) /Z2[0] 

dvdy [0] = (dZdy [1] -Z2[l] *dZdy[0] /Z2[0]) /Z2[0] 

dvdy [1] = (dZdy [2] -Z2 [2] *dZdy[0] /Z2[0]) /Z2[0] 



Viscous Evaluation FOR(cnode,o,3H /* loop on edges of the triangle */ 
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double pia, dGam, Dt; 
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mu_f rom_T ( Temp ) ; 


276 APPENDIX E. NAVIER-STOKES SOLVER AND MESH ADAPTION CODE 






F0R(i,0,3H 

normal [i] [0] = nCoord[cyclic_minusO (i)] [1] - nCoord[cyclic_plusO(i)] [1] ; Form the distributions, 

normal [i] [1] = nCoord[cyclic_plusO(i)] [0] - nCoord[cyclic_minusO(i)] [0] ; 



Finite Volume Functions 
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111 = (double *) calloc ( nodes, sizeof (double) ) 

112 = (double *) calloc ( nodes, sizeof (double)) 
122 = (double *) calloc ( nodes, sizeof (double)) 
lfl = (SolVec *) malloc ( nodes*sizeof (SolVec)) ; 
lf2 = (SolVec *) malloc ( nodes*sizeof (SolVec)) ; 



stack_nodes( nodal_stack ); 
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The limiter is defined as, At that point we can evaluate the limiter, choosing the appropriate one of (U 
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quadrature_point ( node->geom . coord, other_node->geom. coord, 
(celln->cell) ->geom.metric . split_point , QP 
F0R(in,0,NDIMS) r[in] = QP[in] - node->geom. coord[in] ; 



Umin[eq] = UO[eq] + psi[eq] 
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FOR(i ,0 ,NDIMS) r[i] = QP[i] - x0[i] 
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Node *n0, *nl ; /* node pointers */ >* J *' 

SolVec distO, distl ; /* temp distribution vectors */ 

« unsigned int ncheck ; */ > /* end loop on edges of the triangle 



6.7 Finite Volume Artificial Dissipation 
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The projected velocities. V = (w, v) • n. Also, get the speed of sound to left and right, used 
Reconstruct to the left and right states of the conserved variables. Then decode the param- - u eigenvalue limiting 

ctcr, primitive, and Roc-averaged variables at the face. 



veltilde = Vtilde [1] * normal[0] + Vtilde[2] * normal [1] ; X[l] [2] = Vtilde[l] + atilde * normal [0] 
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6.8 Boundary Flux for Finite Volume 

X[0] [0] = X[0] [2] = X[0] [3] = 1 ; 

x[0] [1] = 0 • Evaluate and distribute the flux through a boundary edge. The interior is to the right of 

X[1] [0] = vtilde [l] ; the edge, and the boundary to the left. The quadrature points are at the j and f points 

X[l] [1] = - normal [1] ; along the edge, with distributions going to the nearest node. 



exterior 
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printf ("Bad BC = a / 0 d for edge with end points\n( %f, %f ) , ( %f, %f )\n", case freestream_second: /* second freestream 
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Fluctuation Splitting Functions Node **node = cell->geom. connect. nodes 




dW_dZ( cell->tilde .parameter , dw_dz 
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The full fluctuation is formed as the sum of the '2-D and axi parts, 



288 APPENDIX E. NAVIER-STOKES SOLVER AND MESH ADAPTION CODE 





7.3 Sign of Wavespeed Vector 
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is 1 for 2-D and y for axisymmetric. Z varies linearly over the triangle. Integration 



with, FZ [ j ] [k] [1] -= 0.75 * pibar * df_dz[j] [k] [1] 
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Eigenvalue limiting is added to $ as, 
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CELL_L00P_AT_N0DE (node) -[ - 2 * norm[0] * norm[l] * node->vars .primitive [2] ; 

F0R(iel,0,3) Vf[2] = ( SQUARE (norm [0] ) - SQUARE (norm [1] ) ) * node->vars . primitive [2] 

pia += (( (celln->cell) ->geom. connect .nodes [iel] )->geom. coord [1] ) - 2 * norm[0] * norm[l] * node->vars .primitive [1] ; 

/(( (celln->cell)->geom. connect .nodes [iel] )->bc ? 2 : 4); break; 

++cln ; 



bedge->bc, node->geom. coord [0] , node->geom.coord[l] ); 7.6 FS Eigenvalue Limiting 
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the 2 nodes for this boundary edge */ 



Axisymmetric Source Term Evaluation 
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gensource( cell->geom. metric . area, a, b, distmp ); 
FOR (cnode ,0,3) distrib[2] [cnode] -= distmp [cnode] / 
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with, 



Governing Equations 9.2 Scalar Flux 


296 APPENDIX E. NAVIER-STOKES SOLVER AND MESH ADAPTION CODE 



S is 

>h 43 


S la 3 3 3 3 

c c a 3 c c 

•H *rl *rl *rl "H *H 

*k *t » tt =tfc =tt 


P 

w 

> 


p 

o 

O 


’ S ^ ^ 


l + ii 

^ S. 


° ^ ^ 


J3 


II 

fe. 


d) 

> 

a) 

■p 







F[2] [1] = U[2] * V[2] + V[3] 
F[3] [0] = Z[l] * Z[3] 

F[3] [1] = Z[2] * Z[3] 



df_dz[l] [3] [1] 
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10 H 3 .b 3 .slli Adaption F0R(i t 0,2) (edge->node[i])->adapt. error += edge->error; 


298 APPENDIX E. NAVIER-STOKES SOLVER AND MESH ADAPTION CODE 


§ a> £ a 

| O h' " 

0£ a H <D 

O A <D ho 

I u T3 

I Q) A <U 

a. ho 1 1 

o a) m 

O 0) b0 -H 


A 43 43 A 43 a u 

43 • A 43 ■ O = 43 O 

■ HH ^ 43 • 03 *H = 43 P 

ooi-h - +j u +j 43 • = to 

■ri T) H 3 ffl .H Q, • B S 43 C 

t'd'OPOBiilnid'ri • <0 

+J+J+J IflrH 15 T3 J3 H MU M 

01 m m 5 <h A rt M 43 H 42 -P 

vvvvv= = = = = = = 


fiSBSSSSBBBSC 

•H *rl *H *H *H tH *H -H *H *H *H *H 



void global_error_swap( void ){ node [2] = node_left; 
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testCellLeft . geom. me trie. split _point) ; 
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Update the edge and node errors. 
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EDGE_LOOP_AT_NODE ( node ){ 
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eAC->cell_lef t ; 
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else if( etemp [0] ->cell_left != etemp[2] ->cell_left && 
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Recompute metric terms. 
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Kill the node, a. cell, and two edges. 
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Add the new cells. Hook to the old nodes. Update nodal values that depend upon the cells. nAB->geom.dual_area += cA->geom. metric .area / 3; 



nCD->geom. dual_area += cD->geom. me trie. area / 3; hook_from_node_to_edge ( nCD, eCD ); 
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longest_edge = node->geom.edgecon->edge; int move_node( Node *node , Coord r ){ 
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default : 
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12 Roe Averaging 13 Stonehenge 
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parameter_to_primitive ( Zave, Vtilde ); Node .nptmpO, .nptmpl ; /» temp pointers to nodes »/ 
»Htilde_p = Zave [NEQS-1] / Zave[0] ; Coord DL . /* temporary lengths */ 
*atilde_p = a_from_H_Vel C »Htilde_p, tVtilde[l] >; double dGam ; /* boundary edge length */ 
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If desired, the viscous terms can be computed on a containment-dual mesh instead of the 
median- dual. 



14 Command Line Options strcat ( strcpyC inrstf, rootf ), 
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15 Read Input Deck else if ( !strcmp( value, 
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else if C !strcmp( key, "axi_source_type " ) ){ /* source option */ 

if ( ! strcmp( value, "mass.lumped" ) ) axi_source_type = mass.lumped ; else if( !strcrap( key, "adapt. tar get" ) ) /* adaption target 



adapt_target = S2doub( " adapt _t arget " , value ) ; hist_skip_f ile = (int) S21int ( "hist_skip_f ile" , value 
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CFL_v = CFL_v ? CFL_v : CFL ; /* if no viscous CFL, set to inviscid */ 16 Initial Grid 
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Allocate memory for a node and then read nodal coordinates from a formatted input file. 



success, cells 
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printf ("Incorrectly read °/id connecting nodes for cell %u\nExiting\n" 



17 Starting Line-Up 
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18 Thermodynamic Routines 18.2 e(p, P) 
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19 Transformations zm = ucn / zco] 
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19.5 AW/AZ 19.6 AU/AW 
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void gradZ_to_gradV( SolVec gradZ [NDIMS] , SolVec gradV [NDIMS] , SolVec 
double dV_dZ [NEQS] [NEQS] ; 
double fac = GAMM1 / GAMMA; 



20 Triangle math functions 
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N0DE_L00P_F0RWARD node->geom . dual_area = 0 ; 

CELL_L00P_F0RWARD FOR (n, 0,3) 

(cell->geom. connect .nodes [n] )->geom. dual_area += cell->ge ora .metric . area/3 ; 



the split-point of a triangle for a containment dual 
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F0R(i,0,nq) F0R(j,0,2) double norm_cross_edges ( Edge *edgel , Edge *edge2 M 

for( k=0, grad_u[j] [i] =0 ; k<3 ; ++k ) double denom; 

grad_u[j] [i] += u[k] [i] * n[k][j] ; denom = dnrm2(NDIMS,edgel->r ,1) *dnrm2 (NDIMS ,edge2->r , 1) ; 

return fabs (cross_edges (edgel ,edge2)) /denom; 



21 West World Re-write the grid file, just in case we adapted it. 
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node->vars . wall[l] * nondira_rho * SQUARE (nondim_V) , ifC next_norra < 0 ) 

node->vars . wall[2] * nondim_rho * SQUARE ( no ndim_ V) , Aetna_node ( node->links .node_number , "\n** Negative resid!" 

node->vars.wall[3] * nondim_rho * pow( nondim_V, 3 ), if( next_norm > max.norm ){ 

node->bc ) ; max_norm = next_norm ; 



Aetna_node( max_node->links .node_number , "\nMax resid at this node"); F0R(i,0,4) fprintfC outplot, " °/ 0 f " , node->vars . wall [i] ); 


328 APPENDIX E. NAVIER-STOKES SOLVER AND MESH ADAPTION CODE 


P "O 

u> a 


=• g 


<H P = = 

- P O 

P fi H - » 

O -ft Oh P P 

rH Jh P O O 

Q. ft 3 H rl 

p <h o a- a. 


p o o 


I «J Q -H 
I T3 « ft 


<u J4 m u, m 


ft 03 ft ft 

eft o cl. a. 

fH ft «H *H 


U _J □ 
I-J u. 
P UJ 
W U 


4^ J>> 

IB 3 
„ o 


in 

E 

-Lft 

0 

E 


§1 


3 H -S 
rt H 


rH 


H '3 

9 ‘3 .2 


rt s: -a 

s >> § 

£ -3 a 

OO g* 


• - p 
<u § 


ft ftJ 
3 J ' 
u u. : 


3J 

d 

o 

a 


| 

h0 

Oh 


£ o 


o 

C/1 


H SfJ 

H .§ .5 


O 

a 

H 


o £ 


fM cn 


p to 1 
o U > 

H J = 


O / 
=S 3 
> = 


3 Q ft = 


a) 3 = 3 = 


h H ft h 1 


ft ft ft ft 


UJ O <U fij 


I*. ft ft fi 


ifl hO ho T3 P «J 
P <U <U o S Jh 


> O «H O ' 


I 03 ft 3 03 


□ (ft, O &4 



"tau-y\" \"Qw\" " L2, Cdi, Cdv, Cdi+Cdv, Cli, Civ, Cli+Clv, Qtot 
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if( bedge->bc == wall_viscous ){ 22 LoCcll Ilicllld© FilGS 
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#define BASICS_head 



#include <time.h> /* time_t */ •{ NGeom geom; Vars vars; Be be; Update update; NFVol fv; 
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file for Boundary Conditions 
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void FV_distrib_cell( Cell *cell ); #define T_max_caloric_perf 800 /* calorically perfect limit, 
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latex: $(SRCl).dvi 



hardcopy: latex $(SRCl).dvi : $ (SRC1 : =. tex) $ (SUB1 : = . tex) $(INCl:=.h) 
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[Swe84] P. K. Sweby. High resolution schemes using flux limiters for hyperbolic conser- 
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Appendix F 

Non-linear Conservation Law 
Solver 


The source code, written in Fortran, is presented for the fluctuation splitting and 
DMFDSFV algorithms to solve scalar conservation law problems. Non-linear advec- 
tion and linear advection-diffusion can both be treated in two dimensions. The grid 
pre-processor code follows the solver code. 
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If this is a finite volume computation, need to compute Vu for second- order 
accuracy. Grad u is lagged at the previous time level in the case of a Gaufi- 
Seidel update. 

Grad u could be skipped periodically as the solution nears convergence. 



Advance one time-step by looping over all nodes and computing the update. If if ( irelax .it. l ) then ! GS or Jacobi timestep 
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(144,*) ic, cfluc(ic) 



Fluctuation Splitter The higer-order scheme is obtained by limiting the fluctuations. 
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Further manipulations can be made by defining, 
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sions only. First zero the summations. 
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(in2e(j ,nedge ,2) + 3)/2),’ Uo=’,Uother, write (6,*) ’Umax:’, Umax,’ Umin: ’ , Umin, ’ U:’, u(j) 
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Inputs Limiter 
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Generalized Superbee limUer[3], ilimiter: 



PQ < 0 Unlimited, ilimiter: 
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Define wavespeed A. Form cell-average advection velocity as linear average of 
nodal velocities. 



i,2,2)*dxyz(edge2 ,1) , ic2e (itri,2,2)*dxyz(edge2,2) , 
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in2c( nodegs, ntri, 1 ) ! Current triangle $ dxyz(edgel ,ndim+l) **2 
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Nodal update 
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The artificial dissipation 



The numerical flux. Viscos 
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else if ( idir .eq. 6 ) then 
speedx = 0. 
speedy = 0. 

Smith and Hutton problem. A = (2y(l — x 2 ), — 2x(l — y 2 )) 



Rburg 
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FLS3=main domain getgradu spacenode 
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-Ogrep -c ’Citation .* undefined.’ $(PR0G).log && bibtex $(PR0G) 
-@grep -c ’Rerun to get cross’ $ (PROG). log && latex $(PR0G) 
-@grep -c ’may have changed’ $ (PROG). log && latex $ (PROG) 
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Geometric computations 
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Input boundary pointer information. | write (6,*) 
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, a2 , ’ , * , a3 , ’ , ’ , a4 , ’ , ’ , a5 in VGR1D format. Galled from ingrid. 



subroutine vgridin write (6,*) ’stopping in Vgridin’ 
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Median-dual areas about each node are formed as distribution sums. Triangles 
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Scaled normal to this section of the median dual boundary 



do nodet =1,2 r(ie,l,l,3) 
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(6,*) ’increase format statement for 



(22,103) ! LaTeX file header stuff areamd Median-dual areas about each node, Sud, (node) 
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